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ABSTRACT:  Many  field  research  projects  have  been  conducted  to  study  the  effects  of  natural  foliage  on  the 
propagation  and  attenuation  of  sound.  This  research  takes  natural  foliage  into  a  controlled  laboratory  setting  to  test  its 
low-frequency  acoustic  characteristics.  Absorption  of  low-frequency  components  of  unwanted  noise  is  of  interest  to  the 
Army,  but  has  been  an  unsolved  problem  due  in  part  to  the  cumbersome  and  expensive  testing  facilities  needed  to  study 
long  wavelengths. 

In  this  research,  low-frequency  absorption  and  reflection  coefficients  were  found  reliably  and  consistently.  Due  to  study 
of  the  steady  state  conditions,  the  methods  presented  here  could  constitute  a  more  consistent  method  than  ever  before. 
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1  INTRODUCTION 

A  multitude  of  research  projects  have  been  conducted  to  study  the  effects  of 
natural  foliage  on  the  propagation  and  attenuation  of  sound.  Scores  of  species  have 
been  tested  and  compared.1  Other  variables  such  as  wind,  temperature,  humidity,  and 
elevation  have,  for  the  most  part,  been  out  of  the  control  of  the  researcher.  This  paper 
takes  natural  foliage  into  a  controlled  laboratory  setting  to  test  its  low-frequency  acoustic 
characteristics. 

Three  mechanisms  that  have  commonly  been  studied  and  quantified  are  ground 
reflection,  scattering,  and  absorption  of  sound.  Absorption  of  low-frequency  components 
of  unwanted  noise  is  of  interestto  the  U.S.  Army  Corps  of  Engineers,  but  is  an  unsolved 
problem  due  to  the  cumbersome  and  expensive  testing  facilities  needed  for  long 
wavelengths  and  parameters  that  are  out  of  the  control  of  the  researcher. 

Until  now,  low  frequencies  of  sound  (20-80  Hz)  have  not  been  studied  in 
controlled  laboratory  conditions  with  testing  apparatuses  the  size  of  a  half  wavelength  or 
more.  Until  Burns2  brought  pine  tree  components  into  the  laboratory,  experiments  were 
done  in  the  field  under  nonideal  conditions.  However,  his  testing  apparatus  was  a  box 
much  smaller  than  the  half-wavelengths  of  low  frequency  sound.  Also,  of  the  papers  I 
researched,  only  a  few  address  low  frequencies.  Still,  even  those  papers  take  into 
account  only  1/3-octave  bands  as  low  as  band  number  18  (centered  at  63  Hz)  while  the 
focus  of  study  tends  to  be  on  much  higher  bands. 

The  proceeding  research  describes  the  main  testing  apparatus,  a  concrete 
waveguide  herein  called  "the  chamber,"  in  Chapter  2.  MATLAB  code  written  for  using 
the  chamber  serves  as  an  operator's  manual  for  future  researchers.  All  of  the  code  is 
included  and  well  documented  in  Chapter  7. 
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This  thesis  takes  an  electrical  approach;  a  pressure-voltage  analogy  is  used  to 
compare  the  acoustic  setting  with  an  electrical  one.  The  waveguide  is  analogous  to  an 
electromagnetic  cavity  resonator  (or  electrical  transmission  line  in  one  dimension),  and 
analogous  parameters  are  laid  out  in  a  table  in  order  to  make  solutions  of  the  wave 
equation  in  either  analogy  straightforward. 

The  quality  factor  Q  of  the  chamber  is  studied  extensively  in  Chapter  3.  Varying 
amounts  of  vegetation  are  placed  inside  the  chamber  to  affect  the  Q.  One  goal  of  such 
tests  is  to  find  a  relationship  between  Q  and  the  absorption  coefficient  a  of  vegetative 
material.  This  relationship  is  still  under  investigation  and  could  be  a  subject  of  future 
research. 

Results  for  experiments  with  bales  of  straw  and  a  full  White  Pine  tree  are 
included  in  Chapter  4.  Gravel,  atone  time  outside  of  the  scope  of  this  paper,  is  tested  in 
order  to  compare  foliage  results  with  a  nonvegetative  absorber. 

In  Chapter  5,  the  electrical  analogy  is  used  to  find  the  impedance  and  reflection 
characteristics  of  anything  placed  in  the  chamber.  Once  again,  straw  bales  and  tree 
branches  are  used  as  test  foliage.  MATLAB  code  is  supplied  in  an  appendix  and  can  be 
used  by  future  researchers. 

Chapter  6  shows  results  of  how  a  standing  wave  inside  our  testing  chamber  is 
affected  by  adding  increasing  amounts  of  vegetation  to  dampen  the  system. 

The  chamber  handbook  is  presented  in  Chapter  7.  MATLAB  code  along  with 
extensive  commenting  is  included  so  that  future  researchers  can  copy  and  comprehend 
the  routines  used  herein. 

Finally,  conclusions  are  given  and  a  discussion  of  the  limitations  of  the  results 
obtained  is  presented.  This  work  does  make  use  of  other  researchers' testing  methods, 
however  it  does  not  attempt  recreate  their  work. 


3 


2  TESTING  FACILITY  AND  OPERATION 

The  testing  facility  is  the  main  component  that  makes  this  research  original  and 
unique.  The  concrete  chamber  described  in  Section  2.1  could  not  be  recreated  without 
a  sizeable  budget,  an  abundance  of  space,  and  long-term  planning.  MATLAB  code  was 
written  to  work  with  and  process  the  signals  taken  on  the  recording  device.  This 
recorder  is  described  in  detail  in  Section  2.2.  Section  2.3  describes  the  approach  taken 
in  reconciling  transmission  theory  and  operating  the  facilities. 

2.1  Impedance  Chamber 

The  main  testing  apparatus  for  this  research  is  a  reinforced  concrete  waveguide3 
inside  the  laboratory  at  the  U.S.  Army  Corps  of  Engineers  Construction  Engineering 
Research  Laboratory  (CERL).  The  outside  of  the  chamber  is  shown  in  Figure  2.1. 


Figure  2.1  Reinforced  concrete  impedance  chamber  with  steel  door. 


4 


The  chamber's  design  is  such  that  its  2  m  x  2  m  cross-sectional  dimensions  are 
much  smaller  than  its  length  of  8.8  m.  Thus,  the  first  four  modes  of  oscillation  of  an 
acoustic  field  inside  the  chamber  are  the  modes  (1,0,0),  (2,0,0),  (3,0,0),  and  (4,0,0) 

corresponding  to  roughly  20,  40,  60,  and  80  Hz,  respectively.  The  current  research 
deals  with  one-dimensional  uniform  time-harmonic  plane  waves. 

Minimal  displacement  occurs  in  the  reinforced  concrete  walls  during  internal 

acoustic  reflection  so  coefficients  of  internal  reflection  r  (the  overbar  denotes  that  r  is 
complex)  approach  an  ideal  value  of  1.0  with  no  phase  shift.  A  door  closes  at  the  front 
of  the  chamber,  allowing  items  to  be  placed  inside  for  experiments.  A  loudspeaker 
shown  in  Figure  2.2  is  sealed  in  a  small  cutout  in  the  center  of  the  wall  opposing  the 
door,  herein  called  the  back  of  the  chamber.  This  loudspeaker  is  the  noise  source  for 
normal-incidence  plane  waves  and  can  input  a  pulse,  white  noise,  pink  noise,  a  pure 
tone,  or  any  other  signal  into  the  chamber. 


Figure  2.2  Loudspeaker  installed  in  the  center  of  the  back  of  the  chamber. 
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Microphone  ports  shown  in  Figure  2.3  are  cut  out  of  the  chamber  wall  at  different 
heights  and  locations  along  its  length.  This  research  makes  use  of  the  ports  that  are  1 
m  above  the  chamber  floor  (halfway  up  the  2  m  wall).  When  a  steady  state  pure  tone  is 
input  to  the  chamber,  these  microphones  record  the  pressure  fluctuations  of  a  standing 
pressure  wave. 


Figure  2.3  The  microphone  port  on  the  left  has  a  cable  running  out  of  it,  whereas  the  other 
microphone  ports  are  plugged. 

2.2  Yokogawa  ScopeCorder 

Although  almost  any  data  recorder  could  be  used  in  these  experiments,  the 
Yokogawa  ScopeCorder  DL750,  which  is  available  at  the  testing  facility,  is  important  for 

four  main  reasons.  First,  the  DL750  can  record  at  sampling  rates  up  to  2  x  106 
samples  per  second.  Although  such  high  sampling  rates  and  frequencies  are  not 
necessarily  needed  for  the  low-frequency  scope  of  this  paper,  they  are  beneficial  for 
testing  code  over  a  large  range  of  conditions.  Using  such  a  versatile  machine  allows 
future  researchers  to  make  use  of  the  code  contained  in  this  paper  even  if  their  research 
is  targeted  at  higher  frequency  bands. 
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Second,  16  channels  are  available  at  one  time,  making  time  sync  between  all 
microphones  in  the  chamber  simple  and  reliable.  Third,  the  ScopeCorder  is  flexible  in  its 
methods  of  downloading  data  digitally  for  use  and  processing  by  computer  programs. 
Specifically,  the  ScopeCorder  has  an  internal  hard  drive,  an  IBM  microdrive  that  holds  a 
1  GB  PC  card,  and  an  Ethernet  port  for  input  and  output  of  data  via  a  category-5  cable. 
Fourth,  this  data-recording  device  has  a  large  full-color  LCD  screen  on  which  to  view  the 
channels'  inputs  before,  during,  and  after  the  data  are  taken.  This  is  incredibly  useful  to 
pinpoint  a  bad  microphone,  preamp,  power  supply,  or  channel  setting  thus  saving  time 
when  debugging  faulty  systems. 


2.3  Theoretical  Operational  Rationale 

Taking  an  electrical  approach  to  an  acoustical  problem  simply  means  using 
electromagnetic  theory  and  parameters  in  the  following  wave  equation,4, 5,6,7  shown  here 
for  acoustic  plane  waves  propagating  in  the  z -direction  where  p  is  pressure,  c  is  the 
speed  of  the  wave,  and  t  is  time: 


d2p  _  1  d2p 

~dz2~^~dir 


(2.1) 


Table  2.1  shows  the  parameters  of  the  "pressure-voltage"  analogy  and  their  units 
used  throughout  the  rest  of  this  paper.  There  are  other  valid  analogies  such  as  the 
"pressure-current"  analogy  that  will  not  be  covered  here. 

In  order  to  understand  plane  wave  propagation  in  the  empty  chamber,  MATLAB 
code  based  on  transmission  line  theory  was  written  and  is  included  in  Appendix  A.  The 
MATLAB  function  run chamber .m  allows  the  user  to  input  a  frequency  of  pure  tone  at 
which  they  would  like  to  excite  the  chamber,  initially  at  rest.  It  makes  use  of  concepts 
such  as  impedance,  reflection,  Ohm's  law,  and  superposition  all  while  outputting  three 
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graphs:  a  positive-going  digital  pressure  wave,  a  negative-going  digital  pressure  wave, 
and  a  superimposed  total  digital  pressure  wave. 


Table  2.1  The  acoustic/electric/mechanical  analogy 


Acoustic 

Electrical 

Mechanical 

Parameter 

Units 

Paramter 

Units 

Parameter  Units 

Compliance 

ax'  nr 

Pa  ~  N 

Capacitance 

r 

N  ■  m 

Compliance 
(stiffness _1) 

r 

kg 

Pressure 

Pa  =  — 
m 

Voltage 

ll 

Force 

ii 

K> 

3 

Volume 

3 

m 

Current 

&■[  ^ 
ll 

Particle 

m 

Velocity 

s 

Velocity 

s 

Volume 

Displacement 

3 

m 

Charge 

q 

Distance 

m 

Inertance5  or 
Acoustic 

Mass4 

N-s2  kg 

5  _  4 

m  m 

Inductance 

V-s 

A 

Mass 

kg 

Impedance 

Pa-s 

3 

111 

Impedance 

Q 

Mechanical 

Impedance 

N 

m/ 

/  C 

Resistance 

Pa-s 

3 

m 

Resistance 

Q 

Viscous 

Friction 

/  <3 

N 

m/ 

/  s 

Figure  2.4  shows  a  single  frame  of  the  moving  picture  output  of  runchamber 
after  the  loudspeaker  has  been  turned  on  at  a  constant  pure  tone  of  21.5  Hz  for  0.5  s 
(these  numbers  were  chosen  arbitrarily).  The  top  graph  is  the  positive-going  wave  Pp 

moving  to  the  right,  the  middle  graph  is  the  negative-going  wave  P  moving  to  the  left, 
and  the  bottom  graph  is  the  total  pressure  Ptot  in  the  chamber.  A  distinct  jump  in 
pressure  in  Pn  can  be  seen  at  a  distance  of  2  m  down  the  chamber  axis  due  to  the  onset 
of  a  reflected  wave  with  the  source-produced  P  .  This  jump  travels  back  and  forth  with 


the  traveling  waves. 
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Figure  2.4  A  single  frame  of  the  output  of  runchamber.  The  actual  output  is  a  slower-than- 
real-time  movie  where  the  viewer  can  track  the  movement  and  reflection  of 
pressure  in  the  chamber. 
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3  Q  (QUALITY  FACTOR) 


3.1  Methods  of  Computing  Q 

In  this  section,  we  investigate  several  methods  of  computing  the  Q  of  the 
chamber.  Historically,  Q  designated  the  ratio  of  reactance  to  resistance  of  a  series 
circuit  element.8  It  applies  to  mechanical,  electrical,  and  acoustic  systems.  The  letter 
"Q"was  arbitrarily  chosen  and  eventually  took  on  the  name  "quality  factor." 

3.1.1  Method  1:  the  3  dB-down  method 

The  definition  of  Q  that  describes  a  system's  frequency  selectivity  is 

Q  = - — - -  (3.1) 

Bandwidth  fx-  f2 


where  fr  is  the  resonant  frequency  of  the  system  and  f  and  f2  are  the  -3  dB  points 
(determined  by  half  power)  on  each  side  of  fr.  Maximum  power  delivered  by  a  source 
to  a  system  occurs  at  fr . 


For  example,  the  magnitude  of  impedance 


of  the  resonant  circuit  shown  in 


Figure  3.1  is  shown  in  Figure  3.2  as  a  function  of  frequency.  At  low  frequencies,  voltage 
across  the  resistor  goes  to  zero  because  an  inductor  is  a  short  at  DC. 


i(t)  =  A  cos  (to/) 
v(/)  =  B  cos  (co t  +  0) 
Figure  3.1  A  parallel  resonant  circuit.9 
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At  high  frequencies,  voltage  also  goes  to  zero  because  a  capacitor  is  a  short  as 
/  -» oo.  Maximum  power  is  transferred  (dissipated)  to  the  resistor  when  the  inductor 

and  capacitor  resonate  each  other  at  f..  Bandwidth,  determined  by  the  frequencies 


where  half  peak  power  is  delivered  to  the  resistor,  is  determined  by  points  where 


is 


1/V2  times  the  peak  value  of 


Z 


on  the  frequency  versus  impedance  curve  shown  in 


Figure  3.2. 

Figure  3.2  shows  an  example  of  a  graphical  representation  of  finding  the  Q  from 


the  bandwidth  of  a  system. 


Figure  3.2  The  bandwidth  of  a  parallel  RLC  circuit  whose  resonant  frequency  is  20  H z  and 
where  Q  - 1 . 

Section  3.2  explains  in  detail  how  the  Q  of  an  acoustic  system  can  be  found  using 
pressure  and  volume  velocity  (to  obtain  impedance  as  a  function  of  frequency)  just  as 
easily  as  the  Q  of  an  electrical  circuit  can  be  found  using  voltage  and  current  (to  obtain 
impedance  as  a  function  of  frequency). 

3.1.2  Method  2:  the  energy  method 

Another  way  Q  can  be  defined  for  drivers,  driver/box  combinations, 
electromagnetic  resonators,  acoustic  cavities,  and  other  analogous  systems  is 
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Q  =  2k  x 


(energy  stored  at  resonance) 


Q  =  (o 


(energy  dissipated  at  resonance) 

(overall  reactive  energy  in  the  enclosure) 
(dissipated  power) 


Q  =  2nf 


time  average  of  the  energy  stored  in  the  system 


Q  =  2rc 

Q  =  co 


rate  of  energy  dissipation 

(oscillating  energy) 

(energy  loss  per  cycle  of  oscillation)  ’ 

peak  value  of  energy  stored  in  1  cycle 
dissipated  power 


(3.2) 


10,  n 


(3.3) 


(3.4)11 


(3.5) 


(3.6) 
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These  definitions  are  mathematically  equal  and  directly  related  to  the  other  two  methods 
presented  here. 

3.1.3  Method  3:  the  logarithmic  decrement  method 

A  third  way  Q  can  be  defined  is  via  the  logarithmic  decrement  method.  The 


decrement  8  is  defined  as  the  ratio  of  the  amplitude  of  any  damped  oscillating  peak  to 
the  amplitude  of  the  next  one.13  The  logarithmic  decrement  A  then  is  defined  as 


A  =  ln(c>) 


(3.7) 


and 


0  =  Y-  (3.8) 

A 

When  the  input  of  a  steady  state  system  (such  as  the  one  shown  in  Figure  3.1)  is 
shut  off,  its  response  decays  logarithmically  corresponding  to  e~n ,  where  y  is  an 
attenuation  constant  and  t  is  time.  An  example  of  this  is  shown  in  Figure  3.3. 

The  relationship  between  Q  and  y  can  be  found  from  the  following  procedure. 
Let  A0  and  A /  be  the  amplitudes  of  successive  larger  and  smaller  oscillations, 
respectively.  The  decrement  is  defined  as 
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input  is  completely  cut  off  at  time  t  =  0  s.  The  device's  energy  decays 
logarithmically. 


8  =  ^ 
A 


Let  the  period  of  the  system  be  T , 


T  = 


1 

7 


(3.9) 


(3.10) 


where  /  is  the  frequency  of  the  oscillation.  The  logarithmically  decaying  envelope 
defines  the  amplitudes  of  oscillations  as  follows5 

A  =  An  en  (3.11) 


where  A  is  the  amplitude  of  the  decaying  sinusoid  at  some  time  t.  For  two  successive 
oscillations,  occurring  one  period  T  apart,  Equation  (3.11)  can  be  rewritten  as 


e 


yT 


A 

A 


Substituting  Equation  (3.9), 


0rT  _ 


=  8, 


(3.12) 


Plugging  Equation  (3.12)  into  Equation  (3.7),  we  get 
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A  =  In 


f  r/\ 
e/j 
\  J 


A 


(3.13) 


Finally,  plugging  Equation  (3.13)  into  Equation  (3.8),  we  get 


Q  = 


Y 


(3.14) 


Mathematically,  one  implication  of  this  result  is  that  the  Q  equals  the  number  of  cycles 


undergone  until  the  amplitude  reaches  e  71  (~  0.04)  of  the  initial  value. 

Another  equivalent  perspective  is  to  measure  the  amplitude  of  two  consecutive 

cycles  of  the  decay.  The  second  cycle's  amplitude  will  have  a  value  that  is  e ~*IQ  times 
the  amplitude  of  the  previous  cycle.  If  the  Q  of  the  chamber  is  high,  however,  this 
method  is  difficult  with  empirical  data  because  the  ratio  of  amplitudes  of  successive 
cycles  of  the  decay  will  be  extremely  close  to  1 . 

In  such  instances,  it  is  best  to  count  multiple  oscillations  of  the  system.  If  the 
decay  of  amplitude  in  one  cycle  of  oscillation  were  too  small  to  measure,  one  could 
measure  the  decay  of  some  number  N  cycles  of  decay.  In  this  case,  the  Nth  cycle's 
amplitude  will  have  a  value  that  is  e  /V7r/e  times  the  amplitude  of  the  original  cycle. 


3.2  Choosing  a  Method 

The  first  attempt  to  measure  Q  consisted  of  sending  a  sweepable  pure  tone  into 
the  empty  chamber.  A  microphone  was  placed  inside  the  chamber,  and  its  output  was 
displayed  on  the  screen  of  an  oscilloscope.  This  test  did  not  work  well  because  a  high 
Q  implies  a  very  narrow  bandwidth  (see  Figure  3.2).  The  analog  signal  generator  was 
not  precise  enough  to  dial  in  changes  of  frequency  A/  small  enough  to  find  such  a 
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bandwidth.  The  values  of  the  three  frequencies  needed  for  Equation  (3.1),  namely,  fr, 
fx,  and  f2,  were  determined  by  counting  divisions  on  the  screen  of  the  oscilloscope. 
Although  this  constitutes  a  crude  method,  Q  for  mode  (1,0,0)  was  found  to  be  73  (a 
value  very  close  to  later,  more  precise  tests). 

After  this  initial  testing,  data  was  no  longer  analyzed  in  real  time.  Rather,  it  was 
recorded  and  processed  in  MATLAB  using  signal  processing  routines  to  precisely 
analyze  the  signals. 

The  second  experiment  to  measure  Q  attempted  to  make  good  use  of  the 
aforementioned  logarithmic  decrement  method.  Once  again,  a  signal  generator  with  a 
pure  tone  at  a  user-defined  frequency  was  connected  to  the  chamber's  loudspeaker.  A 
microphone  was  placed  inside  the  chamber  and  this  time  its  output  was  recorded  on  a 
Sony  TCD-D8  DAT  recorder.  The  pure  tone  was  played  long  enough  for  a  steady  state 
standing  pressure  wave  to  be  set  up  inside  the  chamber.  The  input  to  the  chamber  was 
then  abruptly  cut  off,  and  the  decay  of  the  standing  pressure  wave  was  recorded  and 
analyzed.  Q  for  mode  (1,0,0)  was  85 . 

The  third  type  of  test  came  about  in  order  to  measure  more  than  one  frequency  at 
a  time.  For  this  third  group  of  tests,  an  impulse  (or  spark)  was  put  into  the  chamber.  An 
impulse  has  equal  energy  across  all  frequencies.  The  decay  of  the  pressure  wave  due 
to  the  impulse  was  too  difficult  and  erratic  to  analyze  with  logarithmically  decaying 
envelopes. 

The  fourth  and  final  way  to  analyze  the  chamber  was  to  use  white  noise  as  the 
input  source,  once  again  at  steady  state.  White  noise  is  made  up  of  impulses  with 
random  amplitude  in  succession,  so  it  too  has  equal  energy  across  all  frequencies. 
Since  we  are  only  interested  in  the  first  four  resonant  modes  of  the  chamber 
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(frequencies  below  100  Hz),  the  white  noise  was  lowpass  filtered  with  a  cutoff  of  110 
Hz.  Figure  3.4  shows  the  power  spectrum  of  the  filtered  white  noise. 


interest,  i.e.,  15-100  Hz. 

Microphones  along  the  chamber  wall  sent  pressure  signals  due  to  the  white  noise 
inside  the  chamber  to  the  aforementioned  ScopeCorder.  Data  were  sampled  at  500  Hz, 
well  above  Nyquist  frequency  for  signals  in  the  range  of  15-100  Hz.  The  discrete  Fourier 
transforms  (DFTs)  of  the  signals  were  analyzed.  Figure  3.5  shows  the  result  from  a 
single  microphone  after  smoothing  algorithms  detailed  in  Appendix  B. 


10  20  30  40  50  60  70  80 

Frequency  (Hz) 

Figure  3.5  The  magnitude  of  the  DFT  of  a  microphone  in  the  chamber  when  white  noise  is  the 
chamber's  input. 

The  3-dB  down  method  was  used  on  the  peaks  in  Figure  3.5  to  obtain  the  Q  at 


each  resonance.  As  long  as  there  are  enough  points  in  the  DFT  shown  in  Figure  3.5, 
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the  bandwidth  of  system  resonant  peaks  can  easily  be  found.  With  enough  points  so 
that  the  frequency  resolution  Af  =  J/^q  Hz,  this  proved  to  be  a  reliable  and  useful 
method  of  analysis  and  is  the  method  used  throughout  the  rest  of  this  paper. 

3.3  Q  and  a ,  Absorption  Coefficient 

The  absorption  coefficient  a ,  defined14  as  the  ratio  of  the  energy  absorbed  by  a 
surface  to  the  energy  incident  upon  a  surface,  of  an  absorbing  material  placed  against 
the  door  so  as  to  terminate  the  chamber  is  a  dimensionless  constant  related  to  Q. 
Beranek14  equates  a  to  a  measure  of  the  absorbing  power  of  unit  area. 

In  the  past,  measurement  of  a  has  been  a  detailed  process  often  performed  by 
laboratories  created  specifically  for  such  experimental  determination.  The  method  most 
often  used  consisted  of  measuring  reverberation  times  in  a  reverberation  chamber.  The 
research  in  this  thesis  improves  upon  their  limitations  because  it  deals  with  steady  state 
signals,  thus  minimizing  error.  The  particular  relationship  between  a  and  Q  is  currently 
under  investigation  and  may  be  a  topic  of  future  research. 

3.4  Laboratory  Techniques 

Vegetation  placed  inside  will  lower  the  Q  from  that  of  an  empty  chamber, 
however  the  change  is  dependent  upon  the  location  of  the  vegetation.  This  problem  can 
best  be  studied  via  the  electric  analogy. 

In  transmission  line  theory  dealing  with  lossless  lines  varying  sinusoidally  at 
steady  state,  the  location  of  maximum  voltage  occurs  at  a  location  of  minimum  current 
(and  vice  versa).11  A  transmission  line  terminated  by  an  open  circuit  at  one  end  and  a 
short  circuit  at  the  other  end  is  shown  in  Figure  3.6  to  illustrate  one  mode  of  standing 
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waves  and  depict  boundary  conditions  that  go  along  with  such  line  terminations.  Using 
Ohm's  law  (Z  =  v/ 1 ),  an  added  impedance  (no  matter  how  large)  placed  at  a  location 
where  current  is  zero  does  not  affect  the  standing  wave  pattern.  One  could  open  circuit 
the  line  at  current  nulls  and  short  the  line  at  voltage  nulls  without  changing  these  natural 
oscillations.11 


Figure  3.6  Standing  current  and  pressure  waves  shown  on  a  hypothetical  transmission  line 
with  characteristic  impedance  Z0(a).  The  placement  of  even  an  infinitely  large 

impedance  at  a  current  null  does  not  affect  the  standing  waves  (b). 

In  the  acoustic  analogy,  maximum  volume  velocity  in  a  transmission  medium 
occurs  where  the  pressure  is  minimum.  An  impedance  of  any  magnitude  (but 
infinitesimal  thickness)  added  where  the  volume  velocity  is  zero  (maximum  pressure)  will 
not  affect  the  acoustic  standing  waves.  Therefore,  to  study  the  effect  of  added 
impedance  on  the  standing  wave  patterns,  it  is  necessary  to  add  such  impedances  (i.e., 
vegetation  mass)  where  the  pressure  wave  has  a  null.  Due  to  pressure  doubling,  the 
standing  pressure  waves  of  the  empty  chamber  at  any  frequency  will  have  a  maximum 
at  the  chamber  door.  The  shapes  of  the  first  four  fundamental  modes  of  the  empty 
chamber  are  shown  in  Figure  3.7. 

To  best  study  the  effect  of  vegetation  on  the  Q  for  mode  (1,0,0),  the  brush  would 
be  placed  at  location  1/2 ,  where  /  is  the  length  of  the  chamber.  In  the  same  way,  brush 
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would  be  placed  at  1/4 ,  1/2 ,  and  1/8  (or  any  other  nulls  in  the  pressure  waves  shown  in 
Figure  3.7),  for  modes  (2,0,0),  (3,0,0),  and  (4,0,0)  respectively. 


Mode  (1 ,0,0) 


0123456789 

Mode  (3,0,0) 


0123456789 

Mode  (4,0,0) 


0123456789 

Chamber  Axis  (m) 


Figure  3.7  Standing  pressure  waves  in  an  empty  chamber  for  the  first  four  modes.  The  hard 
chamber  ends  are  almost  perfectly  reflecting  and  can  be  modeled  as  an  open 
circuit. 
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4  QUALITY  FACTOR  RESULTS 

4.1  Empty  Chamber 

The  empty  chamber  is  the  canon  to  which  all  other  chamber  conditions  can  be 
compared.  Figure  4.1  shows  the  electrical  analogy  from  which  to  calculate  acoustic 

parameters  in  two  situations.  V s  and  Zs  are  the  AC  source  voltage  and  current, 
respectively. 


Figure  4.1  The  transmission  line  equivalent  of  the  empty  chamber  (a)  and  the  chamber  with 
some  vegetation  placed  at  some  location  along  its  axis  (b).  The  z  -axis  runs  along 
the  centerline  of  the  chamber.  The  voltage  source  represents  the  loudspeaker 
(pressure  source)  in  the  chamber  and  the  open  circuit  at  z  =  0  represents  the 
chamber's  door. 

A  2 -port  network  is  used  to  represent  vegetation  in  Figure  4.1(b).  The  detailed  study  of 
2-port  networks  is  beyond  the  scope  of  this  research,  but  the  methods  of  computing  Q 
presented  in  Chapter  3  do  not  rely  on  such  an  analysis.  Nonetheless,  the  model  is 
presented  here  for  future  research. 

MATLAB  code  attached  in  Appendix  B  makes  use  of  the  pressure-voltage  analogy 
to  obtain  values  for  Q  and  frequencies  for  the  first  four  fundamental  modes  of  plane 
wave  propagation  in  any  test.  Figure  4.2  shows  what  the  microphones  measure. 

The  results  for  the  empty  chamber  are  shown  in  Table  4.1  on  the  following  page. 
These  numbers  are  taken  from  an  average  of  all  of  the  microphones  in  the  chamber. 
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Figure  4.2 


The  microphones  in  the  chamber  measure  pressure  ( Vmic  in  our  analogy)  at 
different  locations.  Only  one  microphone's  measurement  is  shown  here  at  the 
location  zmic.  Theoretically,  the  Q  of  the  chamber  should  be  the  same  for  every 
microphone  along  the  chamber  axis. 


Table  4.1  Empty  chamber  Q 


1  Mode 

Q 

Frequency  (Hz) 

(1,0,0) 

65 

19.64 

(2,0,0) 

114 

39.39 

(3,0,0) 

172 

59.13 

(4,0,0) 

251 

78.80 

Theoretically,  the  Q  for  an  empty  chamber  with  perfectly-reflecting  walls  should 
approach  oo.  Any  absorption  inside  the  chamber  would  yield  a  lower  Q.  Hence,  the 
finite  results  in  Table  4.1  show  that  nonideal  conditions  were  present,  although  all 
Q  >>i .  Therefore,  in  the  investigation  of  absorption,  it  is  assumed  that  the  vegetation  is 
the  only  thing  absorbing  energy  inside  the  chamber  (and  the  chamber  itself  does  not 
dissipate  energy). 


4.2  Vegetation  in  Chamber 

Sundry  vegetation  was  placed  in  the  chamber,  and  quantitative  results  are  included 
in  this  section.  Each  time  something  is  placed  inside  the  chamber,  its  electrical  analogy 

model  can  be  illustrated  as  Figure  4.2,  where  Zv,egeta;/0„  is  some  unknown  impedance  due 


to  the  vegetation. 
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4.2.1  Apple  tree  branches 

A  29 -kg  dry  leafless  branch  of  an  apple  tree  was  placed  in  the  center  of  the 
chamber  in  an  attempt  to  affect  the  Q  of  the  chamber.  The  branches  had  no  effect 
(probably  due  to  the  low  mass  and  lack  of  leaves)  and  the  results  will  not  be  included 
here. 

4.2.2  Bales  of  straw 

Dense  brush  more  robust  than  the  apple  tree  branches  was  needed  to  affect  the 
Q,  so  bales  of  straw  were  used.  Each  bale  of  straw  had  a  mass  of  39.7  kg  and 
measured  0.375  m  tall  x0.5mwide  xlmlong.  Twenty  bales  were  procured  in  order 
to  create  vegetation  walls  of  different  sizes  in  various  locations  along  the  chamber's  axis. 
The  resultant  Q  for  each  resonant  mode  is  a  function  of  two  variables,  namely  the 
location  and  the  amount  of  added  vegetation.  For  this  reason,  tests  with  equally  sized 
straw-bale  walls  were  repeated  with  the  wall  in  different  locations  along  the  chamber 
axis. 

Each  bale  was  situated  so  that  the  0.5  m  width  was  aligned  with  the  z-axis;  thus 

the  area  facing  the  loudspeaker  was  0.375  m  x  lm  =  0.375m2.  The  results  for  straw- 
bale  tests  are  shown  in  Figure  4.3. 

As  predicted  by  the  transmission  line  analogy,  vegetation  placed  at  // 2  affects  the 
Q  for  resonant  modes  (1,0,0 )  and  (3,0,0)  more  than  at  any  other  location  (see  Figure 
4.2(a)  and  (c)).  Likewise,  the  mode  (2,0,0)  Q  is  affected  most  substantially  by 
vegetation  placed  at  // 4  (Figure  4.1(b)).  Figure  4.2(d)  equally  supports  that  mode 
(4,0,0)  Q  is  most  substantially  affected  by  vegetation  placed  at  // 8,  a  standing  wave 


pressure  minimum  for  that  mode. 
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Q  of  the  Chamber  at  Four  Different  Frequencies 


□  Empty  chamber 

j  1  bale  of  straw  in  chamber 
|  4  bales  of  straw  in  chamber 
I  8  bales  of  straw  in  chamber 


*  Vegetation  at  (12  along  chamber  axis 
t  Vegetation  at  f/4  along  chamber  axis 
$  Vegetation  at  (18  along  chamber  axis 


Figure  4.3  For  each  resonant  mode,  the  Q  of  the  chamber  decreases  with  increasing 
vegetation  amounts  placed  inside.  Referencing  the  symbols  at  the  lower  right,  the 
placement  of  the  vegetation  affects  the  modes  differently  as  well  (see  Figure  3.7). 


Of  interest  to  some  readers  and  future  researchers,  Table  4.2  lists  the  numbers  for 


the  graphical  representation  shown  in  Figure  4.3.  Values  where  vegetation  is  placed  to 
theoretically  decrease  Q  the  most  are  shown  in  bold  font. 
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Table  4.2  Chamber  Q  values  with  varying  amounts  of  straw  bales 


Vegetation 

Vegetation  Location 

Amount 

1/2 

1/4 

1/8 

Mode  (1,0,0),  f -19.6 Hz 

None 

64.90 

64.90 

64.90 

1  bale 

61.09 

61.09 

62.09 

4  bales 

33.35 

41.00 

55.16 

8  bales 

8.60 

10.45 

22.78 

Mode  (2,0,0),  f -39.4  Hz 

None 

114.37 

114.37 

114.37 

1  bale 

104.52 

82.95 

99.54 

4  bales 

76.69 

30.43 

42.55 

8  bales 

53.84 

10.03 

8.40 

Mode  (3,0,0),  f  -59.2  Hz 

None 

172.15 

172.15 

172.15 

1  bale 

148.03 

151.31 

141.81 

4  bales 

39.65 

59.83 

35.65 

8  bales 

10.24 

27.12 

11.64 

Mode  (4,0,0),  f -78.8  Hz 

None 

251.49 

251.49 

251.49 

1  bale 

183.88 

195.99 

167.50 

4  bales 

94.62 

103.55 

41.72 

8  bales 

46.99 

64.05 

21.49 

4.2.3  Northern  white  pine  tree 

A  freshly-cut  Pinus  strobes  (White  Pine)  was  brought  into  the  laboratory  for 
testing.  The  specimen  came  from  the  west  side  of  a  15 -m  thick  stand  of  White  Pine  at 
the  Vermillion  River  Observatory  in  Danville,  IL.  It  was  tested  in  the  chamber  within  two 
days  of  being  cut. 

The  9.2 -m  tall  tree  was  cut  down  and  driven  to  CERL  on  the  back  of  a  trailer  in 
two  smaller  sections  totaling  less  than  9.2  m.  The  longest  section  consisted  of  the  top  6 
m  and  130  lb  of  the  tree.  It  was  bushy  with  branches  full  of  needles.  For  this  top  piece, 
the  trunk  measured  5.75  in  in  diameter  at  the  bottom  cut  and  tapered  to  a  single  bunch 
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of  needles  at  the  top.  The  east-facing  side  of  the  tree  had  significantly  fewer  branches 
and  needles  because  it  faced  the  inside  of  the  stand  of  trees  and  received  less  sunlight 
than  the  west-facing  side. 

The  second  section  consisted  of  the  next  1.5  m;  the  115  lb  section  of  trunk  had  a 
medium  amount  of  branches  with  needles.  The  branches  on  the  east  side  were  barren 
and  void  of  needles.  The  trunk  measured  7  in  atthe  bottom  cut. 

Two  bushy  branches  were  cut  from  the  base  of  a  larger  White  Pine  and  piled  on 
the  trailer.  The  branches,  shown  in  Figure  4.4  (a),  weighed  89  lb  and  70  lb, 
respectively.  They  were  added  simply  for  additional  vegetation  mass  in  the  chamber.  In 
addition  to  the  aforementioned  130  lb  and  115  lb  pieces  of  the  tree,  100  lb  of  branches 
that  had  broken  off  during  tree  removal  were  gathered  from  the  tree  and  piled  on  the 
trailer. 


(a)  (b) 

Figure  4.4  Two  branches  (a)  cut  from  another  tree  were  used  for  extra  vegetation  mass  and 
the  top  130  lb  of  the  tree  (b)  atthe  chamber  door. 

Thus,  the  total  mass  from  the  main  tree  was  (130  lb  +  115  lb  +  100  lb)  345  lb, 

with  an  additional  (89  lb  +  70  lb)  159  lb  of  the  two  branches  from  another  tree.  The 


lowest  1.5  m  of  the  main  tree  had  few  branches  and  needles  and  was  left  in  the  forest. 


25 


The  main  130  lb  section  of  the  tree  was  loaded  into  the  chamber  first,  propped 
on  two  sawhorses,  and  is  shown  in  Figure  4.4  (b)  on  the  previous  page.  The  115  lb 
section  of  the  trunk  was  placed  in  the  rear  of  the  chamber  (near  the  speaker). 

Unlike  the  bales  of  straw  in  Section  4.2.2,  it  was  impossible  to  place  the  tree  in 
discreet  locations  of  standing  wave  pressure  minima.  Results  for  four  different  amounts 
of  tree  are  shown  in  Figure  4.5. 


Q  of  the  Chamber  with  Four  Different  Amounts  of  White  Pine 


Mode  (2,0,0),  f~39.4  Hz 


Mode  (4,0,0).  f-78.8  Hz 
(d) 


]  1 43  kg  of  White  Pine  in  chamber 
|  1 63  kg  of  White  Pine  in  chamber 
|  1 88  kg  of  White  Pine  in  chamber 
|  229  kg  of  White  Pine  in  chamber 

Figure  4.5  Indeed,  the  Q  is  a  function  of  the  mass  of  the  vegetation  in  the  chamber. 


For  each  of  the  first  four  modes,  the  higher  amounts  of  tree  in  the  chamber 
resulted  in  lower  Q,  as  expected.  It  is  difficult  to  determine  which  parts  of  the  tree  were 
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absorbing  energy  because  the  tree  did  not  reside  at  a  finite  z  -location.  The  only  way  to 
load  the  tree  into  the  chamber  was  along  the  chamber  axis  where  the  tree  took  up  the 
entire  length  of  the  chamber. 

Once  again,  for  convenience,  the  values  for  Q  are  listed  in  Table  4.3. 


Table  4.3  Chamber  Q  values  with  varying  amounts  of  pine  tree 


Vegetation 

Amount 

(1,0,0) 

Mode 

(2,0,0) 

(3,0,0) 

(4,0,0) 

None 

64.90 

114.37 

172.15 

251.49 

143  kg 

53.62 

94.90 

133.46 

163.69 

163  kg 

49.66 

87.68 

124.82 

160.58 

188  kg 

46.68 

70.78 

119.33 

156.26 

229  kg 

44.99 

64.81 

114.93 

150.72 

4.3  Testing  Limestone  Crushed  Rock 

Absorption  of  low-frequency  energy  is  not  limited  to  vegetation.  The  testing  facility 
and  methods  detailed  here  are  useful  fortesting  other  materials,  such  as  earth  or  gravel, 
as  well.  Three  tons  of  crushed  limestone  rock  were  procured  and  the  container  called 
"the  box" shown  in  Figure  4.6  was  built. 


Figure  4.6  A  2  m  x  2  m  two-level  box. 
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The  lower  level  of  the  square  box  was  built  using  2  x  10s;  hence  it  measures 
9.25"  deep.  The  upper  level  of  the  box  was  built  using  2  x  8s  and  measures  7.25" 
deep.  Each  side  of  the  box  measures  2  m,  equal  to  the  width  of  the  chamber. 

To  continue  with  the  method  of  incrementally  adding  mass  in  the  chamber,  the  first 
step  was  to  place  the  lower  level  of  the  box  inside  the  chamber  as  shown  in  F  igure  4.7. 


Figure  4.7  The  lower  level  of  the  box  placed  atthe  doorside  of  the  chamber. 

A  test  for  Q  was  also  performed  on  the  empty  half-box.  Ideally,  it  would  have 
absorbed  no  energy  so  that  when  the  box  was  filled  with  material  to  be  tested,  the 
decrease  in  Q  could  be  attributed  solely  to  that  material. 

A  test  for  Q  was  performed  on  the  empty  complete  box,  shown  in  Figure  4.8. 


Figure  4.8  The  full  box  placed  atthe  doorside  of  the  chamber. 
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Limestone  mined  in  Kankakee,  IL,  from  the  Vulcan  Materials  Company 
headquartered  in  Lombard,  IL,  was  delivered  to  our  laboratory.  Properties  of  the  CA-7 
gravel  are  shown  in  Table  4.4. 


Table  4.4  Limestone  properties 


1  Property 

Value 

Soundness 

37% 

LA  Abrasion 

28% 

Deleterious  Chert 

0 

Soft  &  Unsound 

0.4 

Clay 

0 

Shale 

0 

Total  Deleterious 

0.4 

Bulk  SPG 

2.658 

SSD  SPG 

2.702 

Absorption 

1.6 

The  rock,  shipped  from  the  Urbana  transfer  yard  in  2004  meets  all  State  of  Illinois 
gradation  and  quality  specifications.  Gradation  information  is  provided  in  Table  4.5. 


Table  4.5  Limestone  gradation 


1  Sieve 

Average 

a(standard  deviation) 

Specification 

1.5 ” 

(37.5) 

100 

0.0 

100-100 

1” 

(25) 

100 

0.1 

90-100 

0.75” 

(19) 

89.2 

2.7 

- 

0.625” 

(16) 

72.4 

4.4 

- 

0.5” 

(12.5) 

48.8 

4.2 

40-56 

0.375” 

(9.5) 

23.8 

4.2 

- 

0.25” 

(6.3) 

7.5 

2.4 

- 

#4 

(4.75) 

4.7 

1.1 

0-10 

#16 

(1.18) 

2.2 

0.4 

- 

#200 

(0.075) 

1.76 

0.38 

0-2.5 

Pan 

(0) 

0.00 

0.00 

- 

The  lower  box,  loaded  with  limestone  as  shown  in  Figure  4.9  was  tested. 
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Figure  4.9  The  lower  half  of  the  box  filled  with  gravel. 

Next  the  full  box,  loaded  with  the  crushed  rock  as  shown  in  Figure  4.10  was  tested. 


Figure  4.10  The 

Finally,  eight  bales  of  straw  were  placed  on  top  of  the  box  full  of  rock  in  order  to 
further  reduce  the  Q.  The  setup  is  shown  in  Figure  4.11. 


Eight  bales  of  straw  on  top  of  the  full  box  filled  with  limestone. 


Figure  4.11 
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Results  for  the  five  box/rock  tests  are  given  in  Table  4.6.  Empty  chamber  Q 
values  are  given  for  comparison. 


Table  4.6  Chamber  Q  values  from  rock  experiments 


Chamber  Condition 

(1,0,0) 

Mode 

(2,0,0) 

(3,0,0) 

(4,0,0) 

Empty 

64.90 

114.37 

172.15 

251.49 

Empty  lower  box 

59.18 

102.84 

169.12 

251.58 

Empty  full  box 

58.21 

99.93 

167.42 

238.44 

Lower  box  (w/gravel) 

57.20 

89.29 

145.02 

171.81 

Full  box  (w/gravel) 

49.39 

74.39 

90.10 

127.10 

Full  box,  gravel,  and  8 
straw  bales 

32.58 

17.83 

27.28 

19.42 

Figure  4.12  shows  the  graphical  representation  of  the  data  in  Table  4.6  for  the 
sake  of  consistency. 
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Q  of  the  Chamber  for  Gravel  Tests 


UlL 

Mode  (2,0.0).  f— 39.4  Hz 
(b) 


180 
160 
140 
120 
100 
80 
60 
40 
20 

Mode  (3.0.0).  f-59.2  Hz 
(c) 

]  Empty  chamber 
|  Empty  lower  box  in  chamber 
|  Empty  full  box  in  chamber 
|  Lower  box  filled  with  gravel 
|  Full  box  filled  with  gravel 
|  Full  box  with  gravel  plus  straw  bales 


Mode  (4,0,0),  f-78.8  Hz 
(d) 


Figure  4.12  The  empty  chamber  compared  to  the  five  other  tests  with  varying  amounts  of  rock 
and  the  box. 
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5  LOAD  IMPEDANCE  AND  REFLECTION  COEFFICIENT 

Transmission  line  theory  allows  us  to  find  the  load  impedance  ZL  and  reflection 

coefficient  T*  of  any  substance  placed  against  the  chamber  door.  Unlike  Figure  4.2, 
Figure  5.1  below  shows  the  electrical  model  with  the  vegetation  terminating  the 
transmission  line. 


z  =  0  z  =  C 


Figure  5.1  The  terminating  load  impedance  is  most  likely  complex. 

The  reflection  coefficient  may  be  complex  (and  introduce  phase  shift  in  the 
reflected  wave)  when  any  surface/substance  takes  the  place  of  the  chamber  door.  The 
amplitudes  of  the  oscillating  pressure  waves  at  each  microphone  can  be  used  to 

determine  standing  wave  ratio  (SWR),  r*(  and  ZL.  The  theory  and  methodology  is 
described  below  and  the  explanation  of  the  MATLAB  code  as  a  "handbook"  to  use  the 
chamber  is  given  in  Section  7.2. 

The  speed  of  sound  is  found  from  the  temperature  Tair,  the  specific  heat  ratio 

(£  =  1.40),  and  the  gas  constant  (Rair  =2.869 xlO2  J/kg/K  at  the  temperature  of  the 

laboratory)  at  that  temperature.6  A  controlled  temperature  of  297  K  was  maintained 
throughout  the  experiments  in  the  laboratory. 

Cair  =  ^IkRaJcnr  =  345 

The  pressure  stated  in  height  of  mercury  mmH  can  be  found  in  the  local  weather 
forecast  or  on  the  internet.  The  specific  weight  of  mercury  yHg  is  133  kN/m3. 


32 


Patm  =rHg-mmHg 

The  static  density  of  air  can  be  found  from  the  following  equation:6 


The  characteristic  acoustic  impedance  inside  the  chamber  can  be  found  and  is  defined 
by:4'7'11 


"7  _  A,  air  ^ air 

s 

u chamber 

where  Schamber  =4  m2  is  the  cross-sectional  area  of  the  chamber. 

In  setting  up  the  signal  generator,  the  user  may  choose  a  frequency  /  at  which 
to  excite  the  chamber  with  a  pure  tone.  The  expression  2nfjcair  is  frequently 

encountered  when  using  the  wave  equation.  It  is  convenient  to  substitute  the  variable 
P  for  said  expression: 


Ohm's  law5,11  says 


U  U,  Ur 


where  P  is  pressure  and  U  is  volume  velocity.  The  subscripts  i  and  r  represent  the 
incident  and  reflected  waves  respectively.  A  simple,  general  solution5,7, 11  to  Equation 
(2.1)  is  the  superposition  of  a  negative-going  (incident)  and  a  positive-going  (reflected) 
pressure  wave  in  the  coordinate  system  presented  in  the  previous  chapter: 

P  =  Pi+Pr 

Letting  the  z  -direction  be  the  direction  of  wave  propagation, 
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P(z,t)  =  Pi  t  +  —  +Pr  t-—  .  (5.1) 

V  air  J  V  ^ air  J 

V  =  U,+Ur  =  P(P-Pr) 

The  expression  corresponding  to  Equation  (5.1)  for  the  phasor  chamber  pressure15  is 

P(z)  =  P,e'j/!z  +  P,e  j/!z 

where  P,  =  P,  ejdp ,  Pr  =  T,Pi,  and  r,  =  TV  . 

7>(z)  =  T](ej/h  +Y,  e  jpz\ 

=  Pi\ej0p  [eiPz  +|f \eidpeJfSz\  (5.2) 

Note  that  the  complex  values  used  in  the  previous  equations  are  not  necessarily  equal  to 
their  similarly  named  counterparts.  For  example,  Pr  is  not  the  same  as  Pr . 

The  real  and  imaginary  parts  of  Equation  (5.2)  are 

Re(P(z))  =  pP  i|cos(0p  +  J3z^  +  Ip,  I  If /•  I  cos +  0P  -  /3z)  (5.3) 

Im(P(z))  =  |p,  sin(6,p  + /3z)+  Pi  T,  sin (dp  +  6r- /3z).  (5.4) 

Because  we  are  interested  in  steady  state  conditions,  we  may  choose  to  let  t  =  0 
when  6P  =0.  All  phase  angles  then  will  be  relative  to  9P.  Setting  0P=  0,  Equations 
(5.3)  and  (5.4)  become 

Re(P(z))=  P  /|cos(/?z)  +  |.P;||rJcos(0r  -  J3z )  (5.5) 

Im(P(z)j  =  Pi  sin(/?z)  +  |p,  sin (0T  -/?z)  .  (5.6) 

Squaring  Equations  (5.5)  and  (5.6): 

Re2  (R(z) j  =  \p,  cos2  (/?z)  +  Pi  r,|  cos2  (Qr  -/?z) 

+  2  Pi  T,  cos(/?z)cos(<9r  -  fiz) 


(5.7) 
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Im2  [P{z)  \  =  Pi  sin^  ( /?z)  +  P,  r,  sin2  (<9r  -  /?z) 

_2_  (5.8) 

+  2  Pi  T,  sin(/?z)sin(#r  -J3z) 

The  magnitude  of  P(z)  can  be  found  by  plugging  Equations  (5.7)  and  (5.8)  into  the 


following: 


P{z)\  =  ^Re2(P(z))  +  Im2(P(z)) 


Re2(P(z))  +  Im2(p(z))  = 

Pi  +  P;  f,  +2  Pi  r,|(cos(/?z)cos(#r -/?z)  +  sin(/?z)sin(#r  -/?z)) 


=  Pi  fl  +  Tr  +  2  Tr  (cos(/?z)cos(#r -/?z)  +  sin(/?z)sin(#r -/?z)n  (5.9) 


An  aside: 


cos ( /?z)  cos ($r  -/?z)  +  sin(/?z)sin(#r  -/?z)  = 
cos(/?z)(cos(#r)cos(/?z)  +  sin(#r)sin(/?z)) 

+  sin(/?z)(sin(#r)cos(/?z)-cos(<9r)sin(/?z)) 

=  cos2  (/?z)cos(<9r)  +  2cos(/7z)  810(6*,-)  sin  (/?z)  -sin2  (/?z)cos(6*r) 


(5.10) 


and  using  the  identity 


cos  u  -  sin"  u  =  cos  2 u 


where  u  =  J3z ,  we  simplify  Equation  (5.10)  to 


cos(/?z)cos(f?r  -/?z)  +  sin(/?z)sin(f?r  -  /?z)  = 
cos (2 flz) cos (0T )  +  2 cos (/?z) sin (6r ) sin (/?z) . 


(5.11) 


Making  use  of  the  identity 


2  cos  (w)  sin  (u)  =  sin(2w) 


where  u  =  /3z ,  we  further  simplify  Equation  (5.11)  to 

cos  (2  (5z)  cos  ( 0V )  +  2  cos  (/?z)  sin  (f?r )  sin  (/?z) 
=  cos(2/?z)cos(#r)  +  sin(2/?z)sin(#r) . 


(5.12) 
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Using  the  identity 


cos  u  cos  v  +  sin  u  sin  v  =  cos 


(m-v) 


where  u  =  2[5z  and  v  =  0T,  we  simplify  Equation  (5,12)  to 

cos (2  /?z)  cos [Gr )  +  sin (2/?z) sin [0r )  =  cos (#r  -  2/?z) . 
Substituting  Equation  (5.13)  into  Equation  (5.9),  we  get 

Re2(p(z))  +  Im2(p(z)) 
i+|r,|  +2|r,-|cos(2/?z-#r)| 

—  I  I  I-  |2 

-Pi¬ 


ts. 13) 


Pi 
P(z) 


1+ 


+  2 


cos (2 J3z-0r) . 


(5.14) 


Equation  (5.14)  is  the  crux  of  the  MATLAB  code  titled  sw_fitter  included  in 


Appendix  C  and  explained  in  Section  7.2.  The  unknown  variables 


Pi 


,  and  0V  are 


determined  by  minimizing  the  error  between  recorded  and  estimated  pressure  levels.16 
To  simplify  the  math,  they  can  also  be  determined  by  minimizing  the  error  E  between 
the  square  of  recorded  levels  p  and  the  square  of  estimated  pressure  levels  defined  by 


j= 1 


N 


e=i :\p/ 

i= i 


Pi 


Pi 


,  2A 


1  + 


1  + 


+  2 


+  2 


cos  ^  2  fiz. .  -Pr) 


cos  ^2  f3z j  —  6?r  ^  j 


(5.15) 


where  pj  is  the  measured  sound  pressure  level  at  the  microphone  port  with  z 
coordinates  z. ,  and  N  is  the  number  of  microphones  (and  thus  channels  on  the 
recorder)  used. 


To  minimize  E ,  let  the  partial  derivatives  of  E  with  respect  to 


Pi 


,  and  6V 


be  zero: 
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This  set  of  nonlinear  simultaneous  equations  can  be  solved  using  the  Newton-Raphson 
method.  Code  used  for  this  research  included  in  Appendix  C  uses  a  least-mean  squares 
(LMS)  algorithm. 

Once  values  for  Tr  ,  0T,  and  Pi  are  found,  P(z)  can  be  calculated  with 

Equation  (5.14).  Next,  SWR  can  be  calculated  with  Equation  (5.19). 4 

maximum  value  of  P(z) 

SWR  = - = - 

minimum  value  of  P(z) 


(5.19) 
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The  variable  r,  can  be  found  by  definition  as  TT  = 


r, 


0j®v 


and  2 


following  equation:4, 11 


Zr 


=  zn 


l  +  r, 

l-fr 


r  can  be  found  by  the 


(5.20) 
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6  STANDING  WAVE  RESULTS 


6.1  Vegetation  Tests 

The  chamber  was  tested  in  a  variety  of  situations  in  order  to  better  understand 
the  effect  of  vegetation.  First,  the  chamber  was  tested  while  empty.  Next,  the  chamber 
was  tested  with  a  3 -bale  wall  of  straw  (bales  stacked  on  top  of  each  other),  a  6 -bale 
wall,  and  a  10-bale  wall  termination  shown  in  Figure  6.1  to  test  the  effect  of  gradually 
increasing  vegetation.  Finally,  the  chamber  was  tested  with  68  lbs.  of  pine  branches 
and  needles  hung  from  the  chamber  door.  This  setup  is  shown  in  Figure  6.1  below. 


(a)  (b) 


Figure  6.1  The  chamber  terminated  with  (a)  10-bale  straw  wall  and  (b)  pine  tree  branches 
hung  on  the  inside  of  the  chamber  door. 

Tests  in  this  section  were  done  at  four  different  pure  tone  frequencies.  The 
frequencies  were  chosen  at  random  within  the  band  of  interest,  namely  20-80  Hz. 
Both  tabular  and  graphical  results  are  possible,  with  an  example  of  the  latter  shown  on 
the  following  page  in  Figure  6.2  for  the  frequency  f4=  70  Hz.  The  curves  are  taken 

from  the  output  of  four  different  runs  of  sw_fitter .m  mentioned  in  the  previous 
chapter,  explained  in  detail  in  Section  7.2,  and  included  in  Appendix  C. 


Pressure  (Pa) 
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-  —  6-bale  straw  wall,  0.5m  thick,  |Fr|  =  0.47,  SWR  =  2.77 

-  10-bale  straw  wall,  0.5m  thick,  |f.|  =  0.16,  SWR  =  1.39 

Figure  6.2  Best  fit  standing  wave  patterns  for  a  pure  tone  input  at  70  Hz  in  a  chamber 
terminated  with  different  amounts  of  vegetation.  Resultant  values  for  terminating 
reflection  coefficient  and  standing  wave  ratio  are  also  given. 


Table  6.1  Reflection  coefficients 


1  Frequency 

Terminating  Condition 

ir,i 

SWR 

h  =  25  Hz 

Empty  chamber 

0.992 

249 

3-bale  straw  wall 

0.990 

199 

6-bale  straw  wall 

0.906 

20.3 

10-bale  straw  wall 

0.914 

22.3 

f2  =  30  Hz 

Empty  chamber 

1.00 

11  200 

3-bale  straw  wall 

0.934 

29.3 

6-bale  straw  wall 

0.932 

28.4 

10-bale  straw  wall 

0.848 

12.2 

h  =  50  Hz 

Empty  chamber 

0.97 

65.7 

3-bale  straw  wall 

0.888 

16.9 

6-bale  straw  wall 

0.740 

6.69 

10-bale  straw  wall 

0.524 

3.2 

U  =  70  Hz 

Empty  chamber 

1.00 

6  720 

3-bale  straw  wall 

0.860 

13.3 

6-bale  straw  wall 

0.470 

2.77 

10-bale  straw  wall 

0.164 

1.39 
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Indeed,  adding  vegetation  decreases  both  the  SWR  and  |r,  |  as  expected.  It  can 

also  be  seen  in  Figure  6.2  that  0T  changes  slightly  for  each  case  by  looking  at  the 
change  in  locations  of  the  peaks  of  the  standing  wave. 

The  case  for  68  lbs  of  pine  tree  branches  (Figure  6.1)  is  not  included  in  the 
following  results  because  it  did  not  differ  from  the  results  of  that  of  the  empty  chamber. 

From  Table  6.1,  it  appears  as  though  straw  absorbs  more  energy  near  the  higher 
frequencies  (relative  to  our  frequencies  of  interest,  20-80  Hz). 

6.2  Gamma  Results 

Another  useful  test  is  to  terminate  the  chamber  with  some  material  of  interest  and 

plot  the  Tr  output  of  sw_fitter .m  as  a  function  of  frequency.  Figure  6.3  shows  the 

results  of  such  a  test  for  three  chamber  conditions. 

The  variance  in  the  top  line  (representing  empty  chamber  conditions)  from 

r,  =1.0  in  Figure  6.3  shows  the  margin  of  error  due  to  noise  in  the  system,  imperfect 

chamber  construction,  and  the  plane  wave  approximation.  The  line  representing  testing 
of  the  pine  tree-branch  wall  shown  in  Figure  6.1  varies  with  frequency  in  an 
unpredictable  manner.  The  wall  of  branches  is  so  nonhomogeneous  that  resonances  and 

reflections  are  unpredictable.  The  smoothest  line  gives  r,  for  the  homogeneous  10- 

bale  straw  wall  shown  in  Figure  6.1.  It  is  shown  that  the  bales  of  straw  absorb  the  most 
energy  in  the  frequency  range  of  70-85  Hz. 
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Magnitude  of  the  Reflection  Coefficient  of  Chamber  Terminating  Conditions 


Figure  6.3 


Tr 


-  68  kg  of  pine  branches  and  needles 
•  10-bale  straw  wall 

as  a  function  of  frequency  for  three  chamber  terminations:  empty,  pine  tree 


sections,  and  bales  of  straw. 
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7  CHAMBER  HANDBOOK 


7.1  White  Noise  Q 

The  code  explained  in  this  section,  titled  white_noise_Q.m  is  used  to  find  the  Q 
of  the  closed  chamber  with  any  vegetation  placed  inside.  It  is  explained  in  a  piecewise 
fashion  here  to  act  as  a  handbook  for  future  researchers  and  is  included  in  its 
uninterrupted  entirety  in  Appendix  B. 

The  header  provides  a  help  menu  for  the  program  and  explains  its  limitations. 
There  is  only  one  input  parameter,  and  that  must  be  a  character  string  with  the  location 
of  a  folder  from  which  the  user  wants  to  process  data. 

function  [frequencies, Qs]  =  white_noise_Q ( folder) 

%WHITE  NOISEQ  Calculates  Qs  of  microphones  in  the  chamber  at  CERL. 

o, 

"o 

%  [frequencies , Qs] =WHITE_NOISE_Q (FOLDER)  makes  use  of  data 
%  taken  on  the  Yokogawa  ScopeCorder  and  prepared  by  code 

%  written  by  Tim  Eggerding  at  CERL.  The  folders  that  hold 

%  the  shots  must  be  renamed  "01",  "02",  etc.  Said  folders 

%  must  be  in  a  main  folder  identified  in  the  character 

%  string  FOLDER.  Valid  for  99  shots  or  less  only. 

g, 

'o 

%  There  can  be  between  1  and  16  channels  of  data  taken  for 

%  any  number  of  shots  (or  events) .  WHITE_NOISE_Q  goes  through 

%  FOLDER  and  averages  the  DFT  of  all  shots. 

g, 

'o 

%  The  averaged  DFT  is  then  smoothed  by  bin  averaging  and  the 

%  Q  is  taken  using  the  "3  dB  down"  bandwidth  method. 

g, 

~o 

%  The  frequencies  at  which  resonant  peaks  are  found  are  also 

%  saved  in  FREQUENCIES. 

g, 

'o 

%  See  also  YOKOGAWADATAGUI . 

%  R.  Lee,  2004 
%  10/01/04 

The  program  loads  pertinent  information  and  has  a  variable  called  bin_size  that  is 
hard-coded  here: 


load  shotTimes  %get  the  time  of  all  of  the  shots 

numshots  =  length (shotTimeArray) ; %total  number  of  shots  to  process 
H_tot  =  0;  %initialize  memory 

bin_size  =22;  %choose  a  bin  size  for  bin  averaging 

%this  parameter  can  be  changed  by  the  user 
%at  a  later  point  in  time  based  on  the 
%sampling  rate  and  length  of  time  of  a 
%single  recorded  shot 


22; 
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A  subroutine  included  in  Appendix  B  titled  bin_average.m  makes  use  ofbin_sizeto 
smooth  out  unwanted  jagged  edges  in  a  waveform. 

If  data  were  taken  at  a  different  sampling  rate  or  for  a  different  length  of  time,  the 
user  would  most  likely  change  this  number  to  something  that  gives  clear,  smooth 
graphical  results.  While  choosing  the  right  number  takes  some  thought,  it  is  best 
determined  by  trial  and  error  and  is  based  on  the  user's  preferences. 

Generating  noise  that  is  perfectly  white  is  difficult.  For  short  recordings,  the 
spectrum  of  the  output  of  the  random  noise  generator  does  not  look  flat  at  all.  In  order  to 
get  noise  that  approaches  true  random  white  noise,  recordings  were  excessively  long, 
increasing  file  size  beyond  what  many  computers  could  process.  Taking  a  larger 
number  of  shorter  recordings  was  the  solution,  and  the  following  for  loop  processes  all 
of  those  recordings  at  once: 

for  k  =  l:num_shots 

%find  out  the  name  of  the  next  folder 

if  k  <  10  %if  it's  a  single-digit  number,  like  "l",  make 
%the  folder  name  a  2-digit  number  like  "01". 
fol_name  =  sprintf ( ' %s%s ' , 1 0 1 ,num2str (k) ) ; 
else  %else  if  it's  between  10-99,  just  name  the  folder 

%the  name  of  that  number, 
fol  name  =  num2str(k); 

end; 

cd ( f olname) ;  %get  into  that  individual  folder 

%get  sampling  rate;  it's  the  same  for  all  channels  on  the  Yokogawa 
cd  CHl;load  eventHeader; cd  . . ;clear  traceName; 
clear  date; clear  time 

%clear  from  memory  parameters  that  aren't  going  to  be  used 

num  channels  =  num  folders (cd) ; %the  #  of  channels  we're  dealing 
%with  =  the  #  of  folders  in  the  working  directory 

%put  together  data  matrix  by  going  into  each  channel ' s  folder 

%and  getting  the  shot 

data  =  build  data (num  channels) ; 

H  =  abs (f f t (data) ) ;  %take  an  DFT  of  the  data  matrix 

d_f  =  samplingRate/size (H, 1) ;  %delta  F,  or  the  smallest  change 

%in  frequency  between  points 
freq_axis  =  0:d  f : samplingRate-d  f; 

%freq_axis  is  not  used  in  future 
%computation,  but  it  is  helpful 
%to  have  it  here  for  plotting 
^purposes,  should  you  want  to 
%view  the  DFT. 

H  av  =  bin  average (H, bin_size) ; 
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save  H  H  %save  variable  for  later  use 

Here  the  discrete  Fourier  transform  (DFT)  of  the  shorter  recording  is  stored: 

save  H_av  H_av 
H_tot  =  H_tot  +  H_av; 
cd  .  . 

end; 

The  DFT,  stored  in  H_tot  ends  up  with  units  that  are  multiple  times  too  large  from  all  of 
the  smaller  shots'  DFTs  being  compounded  on  top  of  one  another.  The  following  line  of 
code  averages  outH_totand  gives  it  proper  units: 

H_tot  =  H_tot/num_shots ;  ^average  'em  out 

%make  a  window 

1  window  =  round (2*samplingRate) ;  %length  in  #  of  samples 

The  variable  i_window  is  a  moving  set  of  samples  that  will  parse  through  the  total  DFT, 
H_tot  to  find  peaks.  Ultimately,  these  peaks  tell  us  the  frequency  of  the  tone  being 
played  in  the  chamber. 

cutoff  freq  =  100;  %cutoff,  in  Hz,  above  which  we 

%don't  care  about  peaks 

A  100  Hz  cutoff  is  instated  so  that  the  code  does  not  have  to  iterate  through  frequencies 
in  which  the  user  is  not  interested. 

num  iterate  =  round (cutoff  freq/df) -lwindow+1;  %#  of  times  our 

%window  will  have  to  be  run  through 
maxx  =  max (H_tot (1 : 1  window, :)); %a  list  of  maxima  from  each  window 

As  the  window  parses  the  DFT,  it  stores  maximum  values  out  of  those  samples  in  the 
window's  range. 

indicator  =  ones (l,num_channels) ;  %initialize  memory 

%has  range  of  [lrlwindow] 

locator  =  ones (1, num  channels) ;  %locates  a  point 

for  k  =  2 : 1  window 

maxx(k, :)  =  max (H_tot (k : k+1  window- 1 ,:)) ; 

%the  maximum  value  amongst  the 
%nearest  samples 

for  imic  =  1 :num_channels 

if  maxx (k, i_mic)  ==  maxx(k-l,i  mic)  %if  the  maximum  in  this 
%set  of  samples  is  the  same  as  the  maximum  of 
%the  last  set  of  samples  (one  window  before) 
indicator (imic)  =  indicator (imic)  +  1; 

%increment  indicator 

else 

indicator (i  mic)  =1;  %if  there's  a  new  maximum, 
end;  %reset  the  indicator  to  1. 


%add  up  the  spectra 
%return  to  parent  directory 
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end; 

end; 

%at  this  point  through  (one  run  through  with  the  window) ,  the  only 
%possible  local  max  could  be  at  zero  Hz  (D.C.) 
for  imic  =  1 :num_channels 

if  indicator (imic)  ==  1  window  %if  so,  | H (1, icolumn) | 

If  a  maximum  value  in  a  window  stays  the  same  as  the  window  incrementally  moves 
across  the  DFT's  samples,  it  must  be  a  local  peak.  If,  however,  new  maxima  occur  as 
the  window  moves,  the  code  is  not  finding  any  peaks.  When  that  occurs,  the  variable 
indicator  is  reset.  The  indicator  is  a  Iso  reset  if  a  local  peak  has  been  determined 
so  that  another  peak  could  be  found: 


local  max (imic)  =  H_tot (1, i  mic) ;  %is  a  local  max 
frequencies (imic)  =  0;  %zero  for  D.C. 

locator (imic)  =  locator (i  mic)  +  l;%increment  it 

Each  time  a  peak  is  found,  locator  stores  the  index  so  thatthe  frequency  of  that  peak 
can  be  determined: 


indicator (imic)  =  1; 

end; 

end; 


%reset  indicator 


for  imic  =  1 :num_channels  %for  each  microphone 

for  k  =  1  window+l:num  iterate  %iterate  through  the  DFT 
maxx(k,i  mic)  =  max (H_tot (k: k+1  window- 1 , imic) ) ; 

%max  value  in  that  window 

if  maxx(k,i  mic)  ==  maxx (k-1, i  mic) %if  this  window's  maxx 
%is  the  same  as  the  last  window's,  increment  the  indicator 
indicator (i  mic)  =  indicator (i  mic)  +  1; 
else  %else  reset  the  indicator 

indicator (i  mic)  =  1; 

end; 

if  indicator (i  mic)  ==  1  window; 

%if  so,  then  | H (1 , i_column) |  is  a  local  max 

When  a  peak  is  found,  it  is  stored: 


local  max (locator (i  mic) , imic)  = 
frequencies (locator (imic) , i  mic) 
%mark  the  frequency  (Hz) 
locator (imic)  =  locator (i  mic)  + 
indicator (i  mic)  =  1; 
indices (locator (imic) -1, i  mic)  = 

end; 

end; 

end; 


H_tot (k, i  mic) ; 

=  f req_axis (k) ; 

1;  %increment  it 

%reset  indicator 
k;  %index  in  the  FFT 


When  peaks  are  known,  then  points  determining  bandwidth  are  found: 


%bandwidth  points  are  at  l/sqrt(2)  times  original  peak 
BW_point  =  local_max/sqrt (2 ) ; 

value  =  localmax;  %a  copy  of  local_max  so  we  can  augment  it 
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With  the  bandwidth  points  found,  the  frequencies  at  those  points  can  be  found.  In  this 
program  they  are  called  f_iowand  f_highto  represent  /2  and  f\ ,  respectively,  of 
Equation  (3.1): 

%get  flow  for  the  "3  dB  down"  method 

for  i  mode  =  1 : size (frequencies, 1) %do  this  for  each  mode  found 
for  i  mic  =  1 : num_channels  %do  it  for  each  channel 

offset  =  0;  %an  offset  (in  samples)  from  f_0 

while  value (i_mode,i  mic)  >  BWpoint (i  mode, i  mic) 
value (i_mode, i_mic)  =  H  tot (indices (i_mode, i_mic) - 

(offset+1) ,i_mic) ; 

offset  =  offset  +  1; 

%increment  the  offset  as  long  as  you're  still  above  the  3  dB  down  point 
end; 

if  indices (imode, i  mic) -offset  >=  1 
%if  the  index  were  ever  recorded, 

%do  some  interpolation  to  calculate  a  low  frequency  to  mark 
%the  bandwidth  point  for  the  "3  dB  down"  method" 
slope  =  Htot (indices (i  mode, i  mic) -of fset+1, i  mic)  - 
H  tot (indices (i_mode, i_mic) -offset, i_mic) ; 
rise  =  BWpoint (imode, i  mic) -H_tot (indices (i  mode, i  mic) - 
offset, i_mic) ; 

f  low (i  mode, i  mic)  =  (rise/slope+indices (i  mode, i  mic) - 
offset) *d_f ; 

end; 

end; 

end; 

%get  fhigh  for  the  "3  dB  down"  method 
value  =  local  max; 

for  imode  =  1 : size ( frequencies , 1) %do  this  for  each  mode  found 
for  i  mic  =  1 : num_channels  %do  it  for  each  channel 

offset  =  0;  %an  offset  (in  samples)  from  f_0 

while  valued  mode,  i  mic)  >  BWpoint  (i  mode,  i  mic) 
value (i_mode, i_mic)  = 

H_tot (indices (i  mode, i  mic) + (offset+1) , i  mic) ; 
offset  =  offset  +  1; 

%increment  the  offset  as  long  as  you're  still  above  the  3  dB  down  point 
end; 

if  indices (i  mode, i  mic) +off set-1  >=  1 
%if  the  index  were  ever  recorded, 

%do  some  interpolation  to  calculate  a  high  frequency  to  mark 
%the  bandwidth  point  for  the  "3  dB  down"  method" 
slope  =  Htot (indices (imode, i  mic) +off set, i  mic)  - 
H  tot (indices (i_mode, i_mic) +of f set-1, i_mic) ; 
rise  =  BWpoint (i  mode, i  mic) - 

H  tot (indices (i_mode, i_mic) +of f set-1, i_mic) ; 
f  high ( imode , i  mic )  = 

(rise/slope+indices (i_mode, i_mic) +of f set-1) *d_f ; 

end; 

end; 

end; 

Equation  (3.1)  is  used  to  calculate  Q\ 

Qs  =  frequencies./ (f  high-flow) ;%the  definition  of  the  3  dB  down  Q 
Return 
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7.2  T  Computation  and  Standing  Wave  Fitter 

The  code  explained  in  this  section  is  used  to  find  information  about  material  placed 
so  as  to  terminate  the  chamber.  The  terms  channel  and  microphone  are  often  used 
interchangeably  since  each  microphone's  output  is  recorded  on  its  own  separate 
channel  of  the  recorder.  Figure  7.1  page  shows  an  example  of  bales  of  straw  that 
were  used  to  test  and  develop  this  code. 

Acoustic  properties  of  interest  are  r, ,  SWR,  and  ZL.  The  code  makes  use  of  a 
pure  tone  put  into  the  chamber  with  a  signal  generator  and  the  microphones  in  the 
chamber  receiving  oscillating  standing  pressure  waves. 

The  frequency  of  the  pure  tone  is  determined,  the  amplitude  of  the  standing  wave 
at  each  microphone  is  calculated,  and  a  standing  pressure  wave  is  "best  fir  to  the  data. 
Figure  7.2  illustrates  one  example  of  the  output  of  the  code. 


T en  bales  of  straw  placed  to  terminate  the  empty  chamber. 


Figure  7.1 
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f  =  46.1  Hz.  |r|  =  0.57,  SWR  =  3.65.  ZL  =  358  -  i  75.7 


O  Measured  microphone  pressure  amplitudes 

Figure  7.2  The  code  outputs  the  best  fit  standing  pressure  wave  as  well  as  the  data  given  at 
the  top  of  the  graphic:  tone  frequency,  reflection  coefficient,  standing  wave  ratio 
(SWR),  and  load  impedance. 

The  code  sw_f  itter  is  explained  in  a  piecewise  fashion.  It  is  included  in  its 
uninterrupted  entirety  in  Appendix  C  for  neatness.  The  program  begins  by  defining 
some  constants: 

%%%  DEFINE  CONSTANTS 

location  rate  =  1000;  %samples  per  meter 

1  chamber  =  8.798;  %length  of  chamber  in  meters 

chamber  axis  =  0 : l/location_rate : 1  chamber; 


This  axis  is  an  x-axis  with  which  to  plot  the  output  shown  in  Figure  7.2: 

T  air  =  297;  %air  temp  in  Kelvin 

R_air  =  2.869E2;  %Gas  constant,  (J/kg/K) 

k  air  =  1.40;  %specific  heat  ratio 

mm  HG  =  762;  %mm  of  HG  for  barometric  pressure;  762  mm  =  30  in. 

gamma_HG  =  133;  %specific  weight  of  mercury  (kN/m^3) 

The  constants  shown  above  are  true  for  laboratory  conditions  at  CERL  in  Champaign, 
IL.  The  following  equation6  is  for  measuring  pressure  relative  to  atmospheric  pressure 
(gage  pressure): 

p  atm  =  gamma_HG*mm_HG;  %atmospheric  pressure 

The  following  equation9  is  known  as  the  ideal  gas  law  and  is  used  to  find  the  density  of 
air: 


rho_air  =  p_atm/ (R_air*T_air) ;  %kg/m^3 
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The  following  equation6  makes  use  of  the  ideal  gas  law  to  find  the  speed  of  sound  in  air, 
cail. ,  at  the  temperature  7\  : 

c_air  =  sqrt (k_air*R_air*T  air) ; %Eq.  1.20  in  Munson 

The  characteristic  acoustic  impedance  is  defined  here: 

Z_0  =  rho_air*c_air/4 ;  %Z  =  rho*c/A  where  A  =  4  m^2 

Loading  the  variable  eventHeader  from  one  of  the  microphone's  folders 
retrieves  the  sampling  rate  of  all  data  taken  for  this  event.  The  Yokogawa  DL750 
ScopeCorder  can  only  record  its  channels  with  one  sampling  rate;  therefore,  retrieving 
the  sampling  rate  for  any  single  channel  suffices  for  all  channels.  In  the  code  below,  the 
data  taken  on  channel  1  in  the  folder  "CHI"  is  used  to  retrieve  the  sampling  rate  for  all 
channels.  If,  for  some  reason,  channel  1  were  not  used  on  the  DL750,  this  code  would 
have  to  be  modified  to  retrieve  the  sampling  rate  from  a  different  channel's  folder. 

%get  sampling  rate 

cd  CHl;load  eventHeader; cd  . . ;clear  traceName; clear  date;clear  time 

The  variables  traceName,  date,  and  time  are  not  needed  for  the  rest  of  the  program 
and  can  be  cleared  to  save  memory. 

loc=  [0.744  1.354  1.964  2.574  3.184  3.794  4.404  5.014  5.624  6.234]'; 
%locations  of  microphones 

The  variable  loc  lists  the  z -coordinate  of  the  first  10  microphones.  It  should  be 
modified  if  a  different  number  of  microphones  are  used  or  some  microphones  are 
omitted: 

%%%  OPEN  THE  FILES  AND  READ  'EM  IN 
num_channels  =  numf olders (cd) ; 

%the  #  of  channels  we're  dealing  with  = 

%  #  of  folders  in  the  working  directory 

The  chamber  has  13  microphone  ports  along  its  length  that  measure 
pressure,  however  Figure  7.2  shows  only  10.  When  all  13  channels  were  used, 
sw_f  itter  gave  an  output  where  the  red  circles  in  Figure  7.2  did  not  necessarily  lie  on 
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the  blue  "best  fit"  standing  wave.  The  error  was  unacceptable  and  can  be  explained  by 
higher-order  evanescent  modes  near  the  loudspeaker. 

Three-dimensional  modes  decay  with  distance  from  the  speaker  as  opposed 
to  the  («,0,0)  modes  (where  n  is  the  number  of  half  wavelengths  shown  in  Figure  3.7) 
that  set  up  standing  waves.  The  math  presented  in  Chapter  5  is  based  on  Equation 
(2.1),  which  describes  a  plane  wave  propagating  in  only  one  direction.  Since  no 
loudspeaker  can  create  a  perfectly  plane  wave,  however,  microphones  closest  to  the 
speaker  pick  up  decaying  higher-order  modes.  Trial  and  error  found  that  data  from 
these  microphones  introduced  the  unacceptable  error  in  sw_fitter.  Data  taken  with 
only  the  first  10  microphones  yielded  results  that  were  acceptable,  predictable,  and 
agreed  with  one-dimensional  wave  theory: 

%put  together  data  matrix  by  going  into  each 
%channel ' s  folder  and  getting  the  shot 
data  =  build_data (num channels) ; 

%outputs  an  array  of  peak-peak  amplitudes  and  frequencies 

The  subroutine  buiid_data  is  included  in  Appendix  B  and  puts  together  a 
matrix  that  contains  all  of  the  microphones'  data  in  units  of  pascals.  Each  channel's  data 
is  placed  in  columns  lined  up  in  numerical  order. 

%%%  FIND  THE  FREQUENCY  OF  THE  EXCITATION 
[amplitude, freq] =f req_f inder (data, samplingRate) ; 

The  subroutine  freq_finder  takes  the  data  matrix  and  outputs  a  list  of  amplitudes 
(one  for  each  microphone)  and  a  single  frequency.  The  frequency  is  that  of  the  pure 
tone  being  put  into  the  chamber.  The  amplitudes  are  the  peak  values,  in  pascals,  of  the 
sinusoidal  oscillating  pressure  from  the  standing  pressure  wave  at  each  microphone. 
These  amplitudes  can  be  seen  as  red  circles  in  Figure  7.2: 

beta  =  2*pi*f req/c_air;  %in  radians 

%%%  MAKE  SOME  INITIAL  GUESSES 
amp  =  sort (amplitude (:),  1  descend’ ) ; 

%tabulate  all  of  the  various  amplitudes 

%make  the  largest  mic’s  amplitude  your  initial  guess 

abs_Pi  =  amp(l);clear  amp; 


51 


The  variable  abs_pi  is  a  good  initial  guess  for  Pi  and  is  taken  from  the  largest  value 


in  the  array  amplitude.  The  variable  d_pi  represents  A 


,  a  small  increment 


between  guesses  of 


Pi 


%an  incremental  amount  to  differentiate  guesses  for  the  magnitude 
%of  the  incident  pressure  wave 
d  Pi  =  abs  Pi/50; 

%initial  array  of  a  bunch  of  guesses  for  | Pi | 

Pi  =  0.2*abs  Pi:d  Pi:1.2*abs  Pi-d  Pi; 


There  is  also  a  variable  for  A0r,  the  increment  between  guesses  of  6r ,  and  A 


the 


increment  between  guesses  of 


%an  incremental  amount  to  differentiate  guesses  of  theta, 

%the  phase  angle  of  gamma 
d  theta  =  pi/50; 

%initial  array  of  a  bunch  of  guesses  for  theta 
theta  =  d_theta:d_theta:2*pi; 

%an  incremental  amount  to  differentiate  guesses  of  magnitude  of 
%gamma,  the  reflection  coefficient 
d  gamma  =  0.002; 

The  user  can  replace  the  very  first  d_gamma  in  the  following  line  with  a  larger  number  by 


the  user  to  make  the  following  for  loops  run  more  quickly.  Doing  so  would  yield  the 


correct  result  only  if  the  resultant 


were  between  the  new  number  and  1 : 


%initial  array  of  a  bunch  of  guesses  for  | gamma | 
gamma  =  d_gamma : d_gamma : 1 ; 


The  next  three  for  loops  create  a  matrix,  p,  representing  P(z ) 


given 


in  Equation 


(5.14). 


The  values  in  p  come  from  combinations  of  guesses  for 


0T,  and 


at 


each  mic  location  zmic  given  in  the  code  by  the  variable  loc: 


%the  following  3  for  loops  generate  a  4 -dimensional  matrix 
%consisting  of  all  of  the  combinations  of  gamma,  theta, 
%mic  location,  and  Pi  amplitude  you  can  get 
for  i_gamma  =  1 : length (gamma) 
for  i  Pi  =  1: length (Pi) 

for  i  theta  =  1 : length (theta) 

p (itheta, :, i  gamma, i  Pi)  =  Pi(i  Pi)* 

sqrt  (1+gamma  (i  gamma)  ^2+2*gamma  (i  jamma) 
*cos (2*beta*loc-theta (itheta) ) ) 1 ; 

end; 


end; 
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end; 

A  matrix  e  is  created  to  measure  the  error  between  p  and  the  actual  measured  pressure 
amplitudes  at  each  mic  location: 

%create  an  error  matrix 

for  igamma  =  1 : length (gamma) 
for  i  Pi  =  1: length (Pi) 

for  itheta  =  1 : length (theta) 

e (i  theta, i  gamma, i  Pi)  =  amplitude  - 
p (i  theta, : , igamma, i  Pi) ; 

end; 

end; 

end; 

The  error,  e  is  squared  to  make  all  values  positive: 

e_2  =  e.^2;  terror  squared 

clear  e; 

A  matrix,  error  is  created  after  summing  values  of  e_2  overall  microphones  (values  of 
zmic).  In  other  words,  each  entry  in  error  contains  the  sum  of  all  of  the  microphone's 


errors  for  that  particular  guess  of 


Gr,  and 


Pi 


It  will  be  one  dimension  less  than  e 


or  e_2  since  those  matrices  have  entries  for  each  individual  microphone: 

error  =  sum(e_2,2);  %create  an  error  matrix  that's  1  dimension 

%less  than  p  summing  10  individual  errors 

clear  e_2 ;  clear  p 


The  minimum  value  of  error  occurs  at  some  particular  guess  of 


0r,  and 


Since  error  does  not  hold  the  actual  value  of 


0r,  or 


Pi 


,  the  following  code  finds 


the  index  (the  location  in  the  matrix  error)  of  the  minimum  value: 

[throw  away, theta  index]  =  min (min (min (error,  []  ,4) ,  []  ,3) )  ;  %throw 
away  the  actual  value,  but  keep  the  index 

[throw  away, Pi  index]  =  min (min (min (error) ,  []  , 3 ) )  ;  %throw  away  the 
actual  value,  but  keep  the  index 

[throwaway, gamma  index]  =  min (min (min (error) , [] ,4) ) ; 

Once  the  index  for  each  of  those  variables  is  found,  the  actual  value  for  the  variable  can 
be  found  by  locating  the  index  in  the  variable  array  theta,  Pi,  or  gamma,  the  list  of 


guesses  for  6r, 


Pi 


and 


respectively: 


theta  =  theta (theta  index); 


%LMS  result  for  theta 
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Pi  =  Pi (Pi  index); 

gamma  =  gamma (gamma  index) ; 


%LMS  result  for  Pi |  amplitude 
%LMS  result  for  gamma | 


Once  the  values  are  found  for  the  variables  that  will  give  the  best  fit  curve  for  P(z ) 


p  is 


recalculated,  this  time  as  a  single  waveform  rather  than  a  matrix  of  guesses: 

p  =  Pi*sqrt (l+gamma^2+2*gamma*cos (2*beta*chamber  axis-theta) ) ; 

The  maximum  and  minimum  values  of  the  curve  p  yield  the  SWR  in  the  following  line  of 
code: 


SWR  =  max (p) /min (p) ;  %standing  wave  ratio 
clear  i; 


The  following  line  of  code  gives  a  single  complex  number  to  represent  Tr  from 


and 


0r 


gamma  i  =  gamma*exp (i*theta) ; %def ine  complex  reflection  coefficient 

The  variable  for  r,-  can,  in  turn,  be  used  to  find  the  terminating  load  impedance,  Z, 
from  Equation  (5.20): 

Z_R  =  Z_0* (1+gamma  i) / (l-gamma_i) ;  %kg/nZ2/s 

The  final  step  is  to  be  able  to  create  a  graphical  output  similar  to  Figure  7.2: 

plot (chamber_axis , p) 
hold;plot (loc, amplitude, ' ro ' ) ; 
xlabel ( 'Chamber  axis  (m) 1 ) ; 
ylabel ( 1  Pressure  (Pa)  1 ) ; 


title (sprintf ( '%s%0.3g%s%0.3g%s%0.3g%s%0.3g%s%0.3g%s' , 

1 f  =  1 , f req,  1  Hz,  Gamma  =  ', gamma,  ',  SWR  =  ',SWR, 

',  Z_L  =  1  ,  real  (Z  R)  ,  1  +  '  ,  imag  (Z  R)  ,  '  i  .  '  )  ) 

Each  time  the  code  is  run,  the  results  are  saved  as  MATLAB  .mat  files  as  well  as  a 
graphical  Adobe  Illustrator  .ai  file: 

save  p  wave  p; 

save  chamber  axis  chamber  axis; 
print  -dill  standing  wave; 
return 
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8  CONCLUSION 

A  new  method  of  testing  vegetation  in  a  controlled  laboratory  environment  has 
been  developed  at  the  U.S.  Army  Corps  of  Engineers  CERL.  Low  frequency  absorption 
and  reflection  coefficients  have  been  found  reliably  and  consistently.  Due  to  study  of  the 
steady  state  conditions,  the  methods  presented  here  could  constitute  a  more  consistent 
method  than  ever  before.  While  there  are  innumerable  species  that  were  not  tested,  the 
procedures  described  in  this  paper  serve  as  a  handbook  of  sorts  for  other  CERL 
employees. 

One  suggested  area  of  improvement  would  be  regarding  the  third  attempt  to 
measure  Q  in  Section  3.2.  Instead  of  analyzing  a  decaying  signal  with  all  decaying 
modes  together,  the  signal  should  be  filtered  for  a  single  mode  and  then  studied. 
Although  the  author  believes  this  method  would  work,  it  defeats  the  benefit  of  using  an 
impulse  as  input  to  the  chamber  because  each  resonant  Q  would  be  studied  separately. 

Another  area  of  improvement  would  be  to  incorporate  the  two-port  network 
presented  the  electrical  model  in  Figure  4.1  into  theoretical  calculations.  This  paper  had 
extremely  reliable  results  when  based  off  of  theory,  and  the  acoustic  properties  of 
vegetation  could  be  more  closely  modeled  if  the  parameters  of  a  two-port  network  were 
known. 

The  largest  leap  of  improvement  would  be  to  find  a  definitive  relationship  between 
the  Q  of  the  chamber  and  the  absorption  coefficient  a  of  some  material.  Through  many 
attempts  at  derivations,  discussions  with  other  professors  and  professionals,  and  a 
number  of  different  mathematical  representations  of  the  chamber,  I  believe  that  such  a 
relationship  is  not  far  away. 


56 


A  suggestion  for  improvement  would  be  to  study  the  impedance  characteristics  of 
the  loudspeaker  and  incorporate  them  into  the  electrical  model.  Once  again,  theoretical 
calculations  would  improve  and  make  the  entire  setup  more  reliable  and  easily 
understood. 


57 


APPENDIX  A 

THEORETICAL  STANDING  WAVE  VIEWER 


The  following  MATLAB  code  was  written  to  obtain  the  graphical  results  shown  in 
Figure  2.4.  It  is  well-commented  and  uses  the  most  basic  transmission  equations.  The 
biggest  benefit  of  this  program  is  to  be  able  to  watch  in  slower-than-real-time  traveling 
waves  interfere  to  create  a  standing  wave.  One  can  see  the  system  go  from  the  onset  of 
an  input  signal  all  the  way  to  a  steady  state  condition. 

function  [p] =runchamber ( f_sound, Rdoor , R  src,T_air) 

% [p]  =  RUNCHAMBER ( f_sound, R  door , R  src, T_air) 

q, 

'o 

%  RUNCHAMBER  simulates  a  loudspeaker  with  a  pure  tone  as  the 
%  source  driving  the  impedance  chamber  we  have  at  CERL. 

o, 

"o 

%  f_sound  is  the  frequency  of  the  sound  source,  in  Hz. 

%  R_door  is  the  reflection  coefficient  at  the  door  (load)  end 

%  R  src  is  the  reclection  coefficient  at  the  source  wall 
%  T_air  is  air  temperature  (use  297  if  none  given)  in  Kelvin. 

o, 

"o 


if  nargin  >  4 

error ( 1  Too  many  inputs .  1  ) 

end; 

if  nargin  <=  3 

Tair  =  297;  %degrees  Kelvin 

end; 

if  nargin  <=  2 


R  src  =  0.95;  %ref lection 

end; 

if  nargin  ==1 

Rdoor  =  0.99;  %ref lection 

end; 

if  nargin  <1 

error ('No  input  arguments! 

end; 

%definition  of  variables 
R  air  =  2.869E2; 
k  air  =  1.40; 

c  jir  =  sqrt  (k  air*R_air*T  air)  ; 
omega  =  2*pi*f_sound; 
beta  =  omega/cair; 
l_chamber  =  8.79; 
dist_rate  =  0.1; 
z  =  0  :  dist_rate : l_chamber ; 

L  =  length(z); 
t_run  =  .5;  % length  of  time 

sampling_rate  =  c_air/dist_rate; 
t  =  l/sampling_rate : l/sampling_rate 

%of  tim 


coefficient  at  source  wall 


coefficient  at  the  door 


Can  1  1 1  compute .  1 ) 


%Gas  constant,  (J/kg/K) 

%specific  heat  ratio 
%Eq.  1.20  in  Munson 
%omega  =  2*pi*f 
%page  415  in  Rao 

%8.79  m  is  the  actual  measurement 
%the  number  of  meters  per  sample 
%an  array  in  . 1  m  steps  for  distance 
%L  =  #  of  samples  along  chamber 
(in  seconds)  to  run  the  simulation 
%number  of  samples  per  second 
:t_run;  %0.25  seconds  is  the  length 
e  you'd  like  to  simulate  the  chamber 


P_src  =  cos (omega . *t) ; 

Pp  =  zeros (length (t) , length (z) ) ; 


%driving  function  (source  pressure) 
%allocate  memory  for  positive- 
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Pn  =  zeros (length (t) , length (z) ) ; 

%load  the  first  sample 
Pp  (1)  =  P_src (1) ; 
n  =  1; 

while  n  <=  t_run*sampling_rate-l 
endpoint  =  Pp(n,L); 
begpoint  =  Pn(n,l); 

Pn(n,L)  =  endpoint*R  door; 


%going  pressure  wave 
%allocate  memory  for  negative- 
%going  pressure  wave 


%clock  pulse 

%for  each  increment  in  time 
%save  endpoint 
%save  beginning  point 
%load  Pn(L) 


Pp(n,l)  =  begpoint*R  src+P_src (n) ;  %load  Pp(l) 


Pp  (n+1 , 2  :  L)  =Pp  (n,  1 :  L  - 1)  ; 

Pn  (n+1 , 1 :  L  - 1)  =Pn  (n,  2  :  L)  ; 
n  =  n  +  1  ; 
if  mod (n, 10)  ==  0 

subplot (3,1,1) ,plot(z,Pp(n, 
axis([0  l_chamber  -3  3]) 
subplot (3,1,2) ,plot(z,Pn(n, 
axis([0  l_chamber  -3  3]) 
subplot (3,1,3) ,plot(z,Pp(n, 
axis([0  l_chamber  -5  5]) 
xlabel ( 1  Distance  (m)  1 ) 
drawnow 

M(n-l)  =  getframe; 

end; 


%shift  Pp(x)  right 
%shift  Pm(x)  left 
%increment  clock 
%every  10  steps,  make  a  plot 
)  ) 

)  ) 

)  +  Pn (n, : ) ) 


end; 


P  tot 


Pp  +  Pn; 


p  =  max (P_tot ( : , 1) ) ; 
return 


%Total  pressure  standing  wave  is 
%the  superposition  of  the  positive- 
%going  wave  and  the  negative- 
%going  wave. 

%maximum  pressure  achieved 
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APPENDIX  B 

CODE  TO  FIND  Q  OF  THE  CHAMBER 


The  following  MATLAB  code  was  written  to  obtain  results  shown  in  Chapter  4. 
Subroutines  are  also  presented  after  the  main  program. 


B.1  White  Noise  Q 


function  [ frequencies , Qs]  =  white_noise_Q ( folder) 

%WHITE  NOISE _Q  Calculates  Qs  of  microphones  in  the  chamber  at  CERL. 

q, 

~o 

%  [ frequencies , Qs] =WHITE  NOISE  Q (FOLDER)  makes  use  of  data 

%  taken  on  the  Yokogawa  ScopeCorder  and  prepared  by  code 

%  written  by  Tim  Eggerding  at  CERL.  The  folders  that  hold 

%  the  shots  must  be  renamed  "01",  "02",  etc.  Said  folders 

%  must  be  in  a  main  folder  identified  in  the  character 

%  string  FOLDER.  Valid  for  99  shots  or  less  only. 

Q, 

'O 

%  There  can  be  between  1  and  16  channels  of  data  taken  for 

%  any  number  of  shots  (or  events) .  WHITENOISE _Q  goes  through 

%  FOLDER  and  averages  the  FFT  of  all  shots. 

Q, 

'O 

%  The  averaged  FFT  is  then  smoothed  by  bin  averaging  and  the 

%  Q  is  taken  using  the  "3  dB  down"  bandwidth  method. 

q, 

'o 

%  The  frequencies  at  which  resonant  peaks  are  found  are  also 

%  saved  in  FREQUENCIES. 


See  also  YOKOGAWADATAGUI . 


%  R.  Lee,  2004 
%  10/01/04 

load  shotTimes  %get  the  time  of  all  of  the  shots 

num_shots  =  length ( shotTimeArray) ;  %total  number  of  shots  to  process 

H_tot  =  0;  %initialize  memory 

bin  size  =  22;  %choose  a  bin  size  for  bin  averaging 

%this  parameter  can  be  changed  by  the  user 
%at  a  later  point  in  time  based  on  the 
%sampling  rate  and  length  of  time  of  a 
%single  recorded  shot 

for  k  =  l:num shots 

%find  out  the  name  of  the  next  folder 

if  k  <  10  %if  it's  a  single-digit  number,  like  "1",  make 
%the  folder  name  a  2-digit  number  like  "01". 
fol  name  =  sprintf ( ' %s%s ' , ' 0 ' , num2str (k) ) ; 
else  %else  if  it's  between  10-99,  just  name  the  folder 

%the  name  of  that  number, 
fol  name  =  num2str(k); 

end; 

cd( fol  name);  %get  into  that  individual  folder 

%get  sampling  rate;  it's  the  same  for  all  channels  on  the  Yokogawa 
cd  CHl;load  eventHeader ; cd  ..;clear  traceName; clear  date;clear  time 
%clear  from  memory  parameters  that  aren't  going  to  be  used 
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num channels  =  num  folders (cd) ;  %the  #  of  channels  we're  dealing 
%with  =  the  #  of  folders  in  the  working  directory 

%put  together  data  matrix  by  going  into  each  channel 1 s  folder 

%and  getting  the  shot 

data  =  build_data (num channels) ; 

H  =  abs ( f f t (data) ) ;  %take  an  FFT  of  the  data  matrix 

d_f  =  samplingRate/size (H, 1) ;  %delta  F,  or  the  smallest  change 

%in  frequency  between  points 
freq  axis  =  0:d  f : samplingRate-d_f ; 

%freq_axis  is  not  used  in  future 
%computation,  but  it  is  helpful 
%to  have  it  here  for  plotting 
%purposes,  should  you  want  to 
%view  the  FFT. 

%save  variable  for  later  use 

%add  up  the  spectra 
%return  to  parent  directory 


%average  ' em  out 


%length  in  #  of  samples 
%cutoff,  in  Hz,  above  which  we 
%don't  care  about  peaks 
num_iterate  =  round ( cuto f f_freq/d_f) -l_window+l;  %#  of  times  our 

%window  will  have  to  be  run  through 

maxx  =  max (H  tot ( 1 : lwindow, : ) ) ;  %a  list  of  maxima  from  each  window 
indicator  =  ones (1, num  channels) ;  %initialize  memory 

%has  range  of  [l:l_window] 

locator  =  ones (1 , num_channels) ;  %locates  a  point 

for  k  =  2 : lwindow 

maxx(k,:)  =  max (H  tot (k : k+l_window- 1 , : ) ) ; 

%the  maximum  value  amongst  the 
%nearest  samples 

for  i_mic  =  1 : num_channels 

if  maxx (k, i_mic)  ==  maxx (k- 1 , i_mic)  %if  the  maximum  in  this 
%set  of  samples  is  the  same  as  the  maximum  of 
%the  last  set  of  samples  (one  window  before) 
indicator (i_mic)  =  indicator (i_mic)  +  l;%increment  indicator 
else 

indicator (i_mic)  =1;  %if  there's  a  new  maximum, 
end;  %reset  the  indicator  to  1. 

end; 

end; 

%at  this  point  through  (one  run  through  with  the  window) ,  the  only 
%possible  local  max  could  be  at  zero  Hz  (D.C.) 
for  i_mic  =  1 : num  channels 

if  indicator (i_mic)  ==  l_window  %if  so,  | H (1, i_column) | 

localmax (i_mic)  =  H  tot ( 1 , i_mic) ;  %is  a  local  max 
frequencies (i_mic)  =  0;  %zero  for  D.C. 

locator (i_mic)  =  locator (i_mic)  +  l;%increment  it 
indicator (i_mic)  =  1;  %reset  indicator 

end; 

end; 

for  i_mic  =  1 : num  channels  %for  each  microphone 

for  k  =  l_window+l : numiterate  %iterate  through  the  FFT 
maxx(k,i  mic)  =  max(H  tot(k:k+l  window- l,i  mic) ) ; 


H  av  =  bin  average (H, binsize) ; 

save  H  H 

save  H  av  H_av 

Htot  =  H  tot  +  H  av; 

cd  .  . 

end; 

H_tot  =  H_tot/num  shots; 

%make  a  window 

lwindow  =  round (2*samplingRate) ; 
cutoff_freq  =  100; 
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%max  value  in  that  window 

if  maxx (k, i_mic)  ==  maxx (k-1, i_mic)  %if  this  window's  maxx  is 

%the  same  as  the  last  window's  then  increment  the  indicator 
indicator (i_mic)  =  indicator (i_mic)  +  1; 
else  %else  reset  the  indicator 


indicator (i_mic)  =  1; 

end; 

if  indicator (i_mic)  ==  l_window; 

%if  so,  then  | H ( 1 , i_column) |  is  a 
localmax (locator ( i_mic) ,i_mic)  = 
frequencies (locator (i_mic) , i_mic) 
%mark  the  frequency  (Hz) 
locator (i_mic)  =  locator (i_mic)  + 
indicator (i_mic)  =  1; 
indices (locator (i_mic) -1, i_mic)  = 

end; 

end; 

end; 


local  max 
Htot (k, i_mic) ; 

=  freq_axis (k) ; 

1;  %increment  it 

%reset  indicator 
k;  %index  in  the  FFT 


%bandwidth  points  are  at  l/sqrt(2)  times  original  peak 
BW_point  =  local  max/sqrt (2) ; 

value  =  localmax;  %a  copy  of  local_max  so  we  can  augment  it 
%get  f_low  for  the  "3  dB  down"  method 

for  i_mode  =  1 : size ( frequencies , 1) %do  this  for  each  resonant  mode  found 
for  i_mic  =  1 : num  channels  %do  it  for  each  channel 

offset  =  0;  %an  offset  (in  samples)  from  f_0 

while  value (imode, i_mic)  >  BW_point (i_mode, i_mic) 
value (i_mode, i_mic)  =  H  tot (indices (i_mode, i_mic) - (of f set+1) , i_mic) ; 
offset  =  offset  +  1; 

%increment  the  offset  as  long  as  you're  still  above  the  3  dB  down  point 
end; 

if  indices (i_mode, i_mic) -of f set  >=  1 
%if  the  index  were  ever  recorded, 

%do  some  interpolation  to  calculate  a  low  frequency  to  mark 
%the  bandwidth  point  for  the  "3  dB  down"  method" 
slope  =  H  tot (indices (i_mode, i_mic) -of f set+1, i_mic)  - 

H  tot (indices (imode, i_mic) -offset, i_mic) ; 
rise  =  BW_point (i_mode, i_mic) -H  tot (indices (i_mode, i_mic) -offset, i_mic) ; 
f_low ( imode , i_mic)  =  (rise/slope+indices (imode, i_mic) -offset) *d_f ; 
end; 

end; 

end; 


%get  f_high  for  the  "3  dB  down"  method 
value  =  localmax; 

for  i_mode  =  1 : size ( frequencies , 1) %do  this  for  each  resonant  mode  found 
for  i_mic  =  1 : num channels  %do  it  for  each  channel 

offset  =  0;  %an  offset  (in  samples)  from  f_0 

while  value (imode, i_mic)  >  BW_point (i_mode, i_mic) 
value (i_mode, i_mic)  =  H  tot (indices (i_mode, i_mic) + (of f set+1) , i_mic) ; 
offset  =  offset  +  1; 

%increment  the  offset  as  long  as  you're  still  above  the  3  dB  down  point 
end; 

if  indices (i_mode, i_mic) +off set-1  >=  1 
%if  the  index  were  ever  recorded, 

%do  some  interpolation  to  calculate  a  high  frequency  to  mark 
%the  bandwidth  point  for  the  "3  dB  down"  method" 
slope  =  H  tot (indices (i_mode, i_mic) +off set, i_mic)  - 

H  tot (indices (imode, i_mic) +of f set-1, i_mic) ; 
rise  =  BW_point (i_mode, i_mic) -H  tot (indices (i_mode, i_mic) +of f set- 

1 , i_mic) ; 

f_high ( i_mode, i_mic)  =  (rise/ slope+ indices (i_mode, i_mic) +offset-l) *d_f ; 
end; 
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end; 

end; 

Qs  =  frequencies ./( f_high  -  f_low) ;  %the  definition  of  the  3  dB  down  Q 
return 


B.2  Number  of  Folders  Subroutine 

function  [numchannels]  =  numfolders ( folder) 
%NUM_CHANNELS  =  NUMFOLDERS (FOLDER) 

%  Numfolders  finds  the  number  of  folders  within  a 
%  folder. 

O, 

'o 

%  'folder'  is  any  string  that  contains  the  location 
%  of  the  main  folder 

O, 

'O 


Ryan  Lee,  2004 


return_path  =  cd; 

cd ( folder) ; 
d  =  dir; 

num  channels  =  0; 
for  i  =  3:size(d,l) 


if  d (i) . isdir  == 
num  channels 

end; 


%save  the  current  path 

%so  that  you  can  return  to  it  later. 

%make  ‘folder’  your  current  directory 

%initialize  #  of  channels  to  zero 

%start  at  i  =  3  because  i  =  1  is  not  for  a  folder, 
%but  rather  the  "."  in  DOS  commands,  and  the  i  =  2 
%is  for  the  DOS  folder,  which  is  definitely 

%not  of  interest  here. 

1  %if,  out  of  all  of  the  files  and  folders  in 
%our  folder,  we  come  across  a  directory  (folder) 

=  numchannels  +  1;  %add  1  to  num  channels  for 

%every  folder  encountered 


end; 

cd (return_path) ;  %return  to  the  original  directory 
return 


B.3  Build  Data  Subroutine 

function  [data]  =  build_data (num  channels) 

% [DATA]  =  BUILDDATA (NUMCHANNELS) 

%  Puts  together  a  matrix  with  data  from  various  channels 
%  of  the  Yokagawa  ScopeCorder  in  columns. 

Q, 

"O 

%Written  for  handling  data  from  the  Yokagawa  ScopeCorder 
%Version  1.0  completed  13  JUL  2004  by: 

%Ryan  Lee 

parent_dir  =  cd; 

%get  the  sensitivities  of  all  the  mics 
sensitivities  =  build_sensitivities (num channels) ; 
d  =  dir; column  =  1; 

for  i  =  1 : num  channels  %for  each  of  the  channels 

cd ( sprintf ( ■ %s\\%s%s 1 , parent_dir , 1 CH 1 , num2str (i) ) ) ; %get  into  folder 
load  eventWavef orm . mat ;  %this  waveform  is  in  volts! 
data (:, column)  =  waveform*sensitivities (column) ; 

%this  column  of  ‘data1  gets  this  folder's  waveform 
%units  will  be  in  pascals 
if  length (waveform)  ~=  size (data, 1) 

error (' Channel  data  is  of  different  lengths!') 

end; 

clear  waveform  %save  memory  by  clearing  variable 
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column  =  column  +  1; 
cd (parentdir) ; 

end; 

return 


B.4  Bin  Average  Subroutine 

function  [Hout] =  bin  average (H, bin_size) 

% [HOUT] =  BINAVERAGE (H, BIN  SIZE) 

%  Bin_average  takes  any  data  vector  or  matrix  and  smoothes  out  the 
%  curve  by  bin  averaging.  The  program  takes  the  nearest  number  of 
%  points  (dictated  by  bin_size  input  variable)  and  averages  them 
%  to  get  a  value  for  that  particular  point. 

q, 

'o 

%  Hence,  size (H  out)  will  equal  size(H). 

O, 

-o 

%  If  H  is  2 -dimensional ,  data  should  be  in  the  columns  of  H,  not  the 
%  rows . 

o, 

*o 

%  @Ryan  Lee 
%  08/13/04 

bin  size  =  2*round (bin  size/2) ;  %make  it  an  even  number 

Hout  =  zeros (size (H) ) ;  %initialize  memory 

Hout (1 :bin_size/2 , : )  =  H (1 :bin  size/2 ,:) ; 

for  i_column  =  l:size(H,2)  %for  each  signal  in  H 

for  i  =  bin_size/2+l : size (H, 1) -bin_size/2 
%RMS  bin  average 

H  out (i , i_column)  = 

sqrt (mean (H ( i -bin  size/2 : i+bin  size/2 , i_column)  . ^2 ) )  ; 

end; 

end; 

return 
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APPENDIX  C 

CODE  TO  FIND  BEST  FIT  STANDING  WAVE 

The  following  MATLAB  code  was  written  to  obtain  results  shown  in  Chapter  6. 


C.1  Standing  Wave  Fitter  Routine 


function  [Pi, gamma, theta, SWR]  =  sw_fitter (filename, locations) 

%SW  FITTER  Calculates  parameters  of  a  transmission  line  terminated 
%  by  a  resistive  load  as  an  analogy  to  a  standing  wave  in  the 

%  chamber  at  US  Army  CERL .  13  microphone  ports  on  the  side  of 

%  the  chamber  must  be  used  unless  the  user  manually  changes 
%  the  variable  "loc"  on  line  _ . 

q, 

'o 

%  The  chamber  can  be  empty,  which  results  in  a  line  terminated 

%  by  a  short  circuit,  or  may  have  a  material  placed  inside 

%  the  chamber  door.  For  current  thesis  research,  the  material 
%  of  choice  is  vegetation  mass. 

o, 

'o 

%  Parameters  output  are  Pi,  the  magnitude  of  an  incident 

%  pressure  wave,  gamma,  the  magnitude  of  the  reflection 

%  coefficient  at  the  chamber  door  due  to  the  added  impedance; 

%  theta,  the  phase  angle  of  the  reflection  coefficient  (gamma 
%  will  be  complex);  and  SWR,  the  standing  wave  ratio. 


[Pi, gamma, theta, SWR] =SW  FITTER (FILENAME, LOCATIONS) 


%%%  DEFINE  CONSTANTS 

location_rate  =  1000;  %samples  per  meter 

l_chamber  =  8.798;  %length  of  chamber  in  meters 

chamber_axis  =  0 : 1/location  rate : l_chamber; 

T_air  =  297;  %air  temp  in  Kelvin 

R  air  =  2.869E2;  %Gas  constant,  (J/kg/K) 

k  air  =  1.40;  %specific  heat  ratio 

mm_HG  =  762;  %mm  of  HG  for  barometric  pressure;  762  mm  =  30  in. 

gamma  HG  =  133;  %specific  weight  of  mercury  (kN/m^3) 

p_atm  =  gamma_HG*mm HG;  %atmospheric  pressure 

rho_air  =  p_atm/ (R  air*T_air) ;  %kg/m^3 

c_air  =  sqrt (k  air*R_air*T_air) ; %Eq.  1.20  in  Munson 

Z _0  =  rho_air*c_air/4 ;  %Z  =  rho*c/A  where  A  =  4  m^2 

gammaR  =  0.99;  %ref lection  coefficient  at  door  of  chamber 

%get  sampling  rate 

cd  CHI; load  eventHeader ; cd  ..; clear  traceName; clear  date; clear  time 

loc= [0.744  1.354  1.964  2.574  3.184  3.794  4.404 

5.014  5.624  6.234  6.844  7.454  8.064]  ■; 

%locations  of  microphones 

%%%  OPEN  THE  FILES  AND  READ  THEM  IN 

num_channels  =  num  folders (cd) ;  %the  #  of  channels  we're  dealing  with  = 

%  #  of  folders  in  the  working  directory 

%put  together  data  matrix  by  going  into  each 
^channel 1 s  folder  and  getting  the  shot 
data  =  build_data (num_channels) ; 

%%%  FIND  THE  FREQUENCY  OF  THE  EXCITATION 
[amplitude, f req] =f req_f inder (data, samplingRate) ; 
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%outputs  an  array  of  peak-peak  amplitudes  and  frequencies 
timeaxis  =  1/ samplingRate : 1/samplingRate : size (data, 1) / samplingRate; 
beta  =  2*pi*freq/c_air;  %in  radians 

%%%  MAKE  SOME  INITIAL  GUESSES 

%tabulate  all  of  the  various  amplitudes 
amp  =  sort (amplitude (:),  1  descend  1 ) ; 

%make  the  largest  mic 1 s  amplitude  your  initial  guess 
abs_Pi  =  amp(l);clear  amp; 

%an  incremental  amount  to  differentiate  guesses  of  theta 
d_theta  =  pi/50; 

%initial  array  of  a  bunch  of  guesses  for  theta 
theta  =  d_theta : d_theta : 2*pi ; 

%an  incremental  amount  to  differentiate  guesses  of  Pi 
d_Pi  =  abs_Pi/50; 

%initial  array  of  a  bunch  of  guesses  for  Pi 
Pi  =  0 . 5*abs_Pi : d_Pi : (3/2 ) *abs_Pi -d_Pi ; 

%an  incremental  amount  to  differentiate  guesses  of  gamma 
d_gamma  =  0.02; 

%initial  array  of  a  bunch  of  guesses  for  magnitude  of  gamma 
gamma  =  d_gamma : d_gamma : 1 ; 

%the  following  3  for  loops  generate  a  4 -dimensional  matrix 
%consisting  of  all  of  the  combinations  of  gamma,  theta, 

%mic  location,  and  Pi  amplitude  you  can  get 
for  i_gamma  =  1 : length (gamma) 
for  i_Pi  =  1: length (Pi) 

for  i_theta  =  1 : length ( theta) 

p (i_theta, : , i_gamma, i_Pi)  = 

Pi (i_Pi) *sqrt (1+gamma (i_gamma) ^2+2*gamma (i_gamma) *cos ( 
2*beta*loc+theta (i_theta) ) ) 1 ; 

end; 

end; 

end; 

%Create  an  error  matrix  that's  1  dimension  less  than  the  dimension  of  p 
for  i_gamma  =  1 : length (gamma) 
for  i_Pi  =  1: length (Pi) 

for  i_theta  =  1 : length ( theta) 
for  i_loc  =  1 : length ( loc) 

e (i  theta, iloc, igamma, i_Pi)  =  (amplitude (i  loc) - 
p (i_theta, i_loc, i_gamma, i_Pi) ) ^2 ; 

end; 

end; 

end; 

end; 

error  =  sum (e, 2 ); clear  e;  clear  p  %summing  13  individual  errors 
%throw  away  the  actual  value,  but  keep  the  index 
[throw_away , theta_index]  =  min (min (min (error, [] , 4) , [] , 3) ) ; 

%throw  away  the  actual  value,  but  keep  the  index 
[throw_away , Pi_index]  =  min (min (min (error) , [] , 3) ) ; 

%throw  away  the  actual  value,  but  keep  the  index 
[throw_away , gamma_index]  =  min (min (min (error) , [] , 4) ) ; 

theta  =  theta ( theta_index) ;  %LMS  result  for  theta 

Pi  =  Pi (Pi_index) ;  %LMS  result  for  Pi [  amplitude 

gamma  =  gamma (gamma  index);  %LMS  result  for  gamma 

p  =  Pi*sqrt (l+gammax2+2*gamma*cos (2*beta*chamber_axis+theta) ) ; 

SWR  =  max (p) /min (p) ;  %standing  wave  ratio 

%define  complex  reflection  coefficient 
clear  i;gamma _i  =  gamma*exp (i*theta) ; 

Z  R  =  rho_air*c  air* (1+gamma  i) / (1-gamma  i) ; 
plot (chamber_axis , Pi*sqrt ( l+gamma^2+2*gamma* 


kg/m^2/ s 
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cos (2*beta*chamber  axis+theta) ) ) 
hold ; plot (loc, amplitude,  1 ro 1 ) ;xlabel ( 1  Chamber  axis  (m)  ' ) 
title (sprintf ( 1 %s%0 . 3g%s%0 . 3g%s%0 . 3g%s%0 . 3g%s%0 . 3g%s 1 , 1 f  = 

freq, 1  Hz,  Gamma  =  gamma, SWR  =  ',SWR, 

Z_L  =  1 , real (Z  R) , 1  +  ' , imag (Z  R) , ' i . ' ) ) 

return 


C.2  Frequency  Finder  Subroutine 

function  [amplitude, frequency] =freq_f inder (data, sampling_rate) 
%FREQ_FINDER  finds  the  main  frequencies  of  waveforms 
%  and  also  the  amplitudes  of  those  signals  for 
%  that  particular  frequency.  The  signal  can  be  a  cal 
%  signal  or  any  other  signal  with  a  main  resonance  that 

%  can  clearly  be  seen  in  the  form  of  a  spike 

%  in  the  FFT  of  the  data  that  is  in  'data'. 

%  That  main  resonance  must  be  above  16  Hz. 

Q, 

"O 

%  The  data  must  be  set  up  in  columns  (as  opposed 

%  to  rows)  in  the  data  input  matrix,  and  there 

%  can  be  more  than  1  waveform. 

O, 

'O 

%  [amplitude, omega] =freq_f inder (data, sampling_rate) 

q, 

“o 

for  i_column  =  1 : size (data, 2)  %for  each  waveform  remove  DC  offset 
data ( : , i_column)  =  data (:, i_column) -mean (data (:, i_column) ) ; 

end; 

fft_length  =  2^ (nextpow2 (size (data, 1) ) -1) ;  %this  will  be  the  length 
%of  the  FFT  and  the  "-l"  insures  that  our  data  will  be 
%truncated  (rather  than  zero  padded)  in  the  FFT  algorithm 
freq_res  =  sampling_rate/f f t_length;  %frequency  resolution  (Hz) 

%  FIND  THE  FREQUENCY  OF  THE  SIGNALS 

num_sigs  =  size (data, 2) ;  %the  #  of  signals  we're  dealing  with 
H  =  f f t (data, f f t_length) ;  %takes  the  FFT  of  the  columns 

%  IMPORTANT  NOTE;  H(l)  is  really  the  element  of  the 
%  FFT  that  is  at  a  frequency  of  0.  Matlab  does  not 
%  index  any  matricies  beginning  with  0,  it  always 
%  begins  with  an  index  of  1.  Therefore  H(0)  is 
%  improper  in  Matlab. 

%  freq_axis  =  0 : f req_res : sampling_rate- f req_res; 
%the  location  (an  index)  of  16  Hz  on  the  frequency  axis 
pointer  =  (round (16/f req_res) +1) ‘ones (1 , num_sigs) ; 
if  pointer (1)  <  fft_length 

H (1 rpointer (1) , ; )  =  0;  %zero  out  everything  below  16  Hz 

%zero  out  everything  above  -16  Hz 
H ( f f t_length+l-pointer ( 1) : f f t_length, : )  =  0; 

else 

error ('FFT  is  not  long  enough  or  Sampling  rate  too  small.'); 
error ( ' . . .Or  actually  maybe  your 

data  was  in  rows  instead  of  columns.'); 

end; 


for  i_sig  =  l:num_sigs  %for  each  waveform 

%keep  going  as  long  as  max  hasn't  been  reached 

while  abs (H (pointer (i_sig) , i_sig) )  <  max (abs (H ( : , i_sig) ) ) 

%increment  pointer  along  the  freq.  axis  of  the  FFT (data) 
pointer (i_sig)  =  pointer (i_sig)  +  1; 

end; 

%in  Hz,  get  frequency  of  that  point  that's  max  in  the  FFT 
freq(i_sig)  =  pointer (i_sig) * freq_res - freq_res ; 
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end; 

frequency  =  mode(freq);  %finds  the  mode  of  these  frequencies 

%and  eliminates  outliers 

%some  cutoff  frequencies 
%signal  freq  -  a  semitone 

cutoff _1  =  round ( frequency*2^ (- 1/12 ) /freq_res) +1 ; 

%signal  freq  +  a  semitone 

cutoff_2  =  round ( frequency*2^ ( 1/12 ) /freq_res) +1; 

cutoff_3  =  size(H,l)+2  -  cutoff_2;  %the  negative  frequency  of  freq_2 
cutoff_4  =  size(H,l)+2  -  cutof fl;  %the  negative  frequency  of  freq_l 
block  size  =  round ( sampling_rate/f requency) ; %in  units  of  samples/cycle 

%an  easy  way  to  zero  out  all  unwanted  low  and  high  frequencies  from 
%zero  to  one  semitone  below  our  calibration  frequency 
H ( 1 : cutoff  1 )  =  0; 

%from  one  semitone  above  our  calibration  freq.  to  one  semitone  below 
H (cutof f_2 : cutof f_3 , : )  =  0;  %the  negative  cal  signal  frequency 

%from  one  semitone  above  the  negative  cal  signal  frequency  to 
H (cutof f_4 : length (H) ,: )  =  0;  %the  end  of  the  FFT  zero  it  out. 

idata  =  real (if f t (H, size (H, 1) ) ) ;  %a  new  time-based  signal.  This 
%should  be  similar  to  the  original  signal  minus  any  unwanted  noise. 
%Its  mean  should  be  zero 

for  i_sig  =  l:num_sigs  %for  each  waveform 
%  FIND  AMPLITUDES  OF  THE  SIGNALS 
n=l ; 

%block  data  up  and  find  max  and  mins 

for  numblock  =  1 : floor ( fft_length/block  size)  %for  each  block 
pos_pt (numblock)  =  max ( idata ( (numblock- 

1) *block  size+1 : num_block*block  size , i_sig) ) ; 
neg_pt (numblock)  =  min ( idata ( (numblock- 

1) *block  size+1 : num_block*block  size, i_sig) ) ; 

end; 

amplitude (i_sig)  =  (mean (pos_pt) -mean (neg_pt) ) /2 ; 

end; 

return 
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